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1.  Introduction 


The  Rachele-Kilmer  microphysics  model  for  very  low  stratus  clouds  and 
subcloud  regions  is  a  largely  theoretically  based  model,  in  contrast  to  current 
empirically  based  models.  The  Rachele-Kilmer  model  is  described  in  a  detailed 
technical  report,  [1]  which  may  be  consulted  for  assumptions  that  compose  the 
model  and  for  mathematical  details.  Some  profiles  simulated  for  very  low 
stratus  clouds  using  this  model  were  compared  with  profiles  calculated  using 
actual  field  data.  [1] 

Use  of  the  Rachele-Kilmer  microphysics  model  for  very  low  stratus  clouds  and 
subcloud  regions  is  enabled  by  computer  program  MACACASM,  the  source 
code  for  which  has  been  stored  in  file  ~  nkilmer/sources/hp_macacasm.f  for 
use  with  Hewlett-Packard  9000  computers  at  U.S.  Army  Research  Laboratory 
(ARL)  West.  MACACASM  is  an  acronym  for  Microphysics  and  Cluster 
Ascent  Cloud  and  Subcloud  Model.  This  program  and  the  other  programs 
described  in  this  guide  are  written  in  Fortran  77. 

Certain  output  files  may  be  used  with  Mie  calculations  or  stored  Mie  efficiency 
factors  to  simulate  vertical  profiles  of  extinction,  backscatter,  scattering,  and 
absorption  coefficients.  This  capability  requires  use  of  one  or  two  other 
programs  also.  Details  are  given  in  sections  5.0,  5.1,  and  5.2. 

This  guide  for  using  the  computerized  model  is  in  two  parts  to  separate  the 
material  required  by  two  different  types  of  users.  The  first  part 
(sections  2  through  6)  is  intended  for  the  scientist-user.  It  assumes  an 
understanding  of  meteorological  theory  and  experience  with  analysis  and 
modeling  of  meteorological  data,  but  no  programming  skills  beyond  basic 
computer  literacy.  The  second  part  (section  7)  has  been  prepared  for  the 
computer  programmer  who  may  be  developing  other  programs  related  to  this 
model. 

NOTE:  MACACASM  is  a  research  mode  program.  It  is  possible  for  the  user 
to  supply  sets  of  input  values  that  are  not  self-consistent.  It  is  up  to  the  user 
to  provide  logically  consistent  input. 
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2.  Input  Required  to  DeHne  Model 


MACACASM  is  designed  for  interactive  use  and  prompts  the  user  for  requir-^d 
input.  The  program  can  be  run  in  batch  mode  using  a  redirected  input  file  to 
replace  interactive  input.  (See  2.1.4  for  an  example  of  such  a  file.) 

When  using  MACACASM  in  interactive  mode,  the  user  must  indicate  choices 
when  prompted.  The  user  may  specify  additional  input  (without  being 
prompted)  by  supplying  it  in  input  files  .haze. defaults  and  zbase_b_c,in.  If 
these  files  and  the  information  expected  to  be  in  them  are  not  found,  default 
values  will  be  used. 

User  input,  required  to  execute  the  current  models,  includes  the  following; 

•  a  method  of  specifying  concentrations  of  all  sizes  of  dropLts  at  the 
reference  height  and  the  radii  of  the  corresponding  dry  nuclei  (see  below) 

•  specification  of  how  the  droplet  number  concentration  varies  with  height 
(inversely  proportional  to  the  volume  of  the  cluster  sphere  in  the  most 
common  usage) 

•  reference  height 

•  relative  humidity,  pressure,  and  temperature  in  the  cluster  sphere,  and 
temperature  in  the  ambient  environment  at  the  reference  height 

•  type  of  nucleus  (preferably  one  of  the  ten  reported  by  Hanel  and 
Lehmann  [2]) 

•  ascent  rate  and  radius  of  the  cluster  sphere  at  the  reference  height 

•  an  entrainment  parameter 

•  the  height  of  the  top  of  the  cloud. 


9 


In  addition,  the  user  specifies  a  value  that  affects  precision  of  calculations:  the 
height  step  used  in  numerical  integration. 

Optional  input  includes  some  constants  in  the  temperature  lapse  rate  equation 
for  modifying  modeled  behavior  near  the  base  of  an  inversion,  which  is 
typically  modeled  as  occurring  at  the  same  height  as  the  top  of  the  cloud.  If 
the  user  does  not  provide  values  for  those  constants,  default  values  are  used. 

Also  optional  is  the  inclusion  of  a  term  that  is  linear  with  height  above  the 
reference  height  in  the  temperature  lapse  rate  equation  for  the  ascending  cluster 
sphere.  The  term  may  be  used  to  model  effects  not  otherwise  modeled, 
including  radiation  effects.  Positive  values  of  the  constant  coefficient  (alapslin) 
in  the  term  tend  to  decrease  values  for  liquid  water  content  compared  with 
what  they  would  be  at  the  same  height  when  the  term  is  not  present  (or  when 
its  constant  coefficient  is  zero).  The  latter-mentioned  case  has  been  called 
quasi-adiabatic.  Negative  values  of  the  constant  coefficient  tend  to  increase 
corresponding  values  for  liquid  water  content.  If  the  constant  coefficient  is  not 
provided,  it  is  considered  to  be  zero.  Use  of  the  constant  coefficient  is 
discussed  in  a  later  paragraph. 

The  entrainment  parameter  is  used  in  modeling  horizontal  entrainment  from  the 
ambient  environment  at  the  same  height  as  the  cluster  sphere.  If  a  value  of 
zero  is  specified,  there  is  an  option  for  providing  values  of  constants  to  be  used 
in  modeling  turbulent  transfer  of  heat,  momentum,  and  water  vapor.  If  the 
latter  constants  are  not  provided,  they  are  considered  to  be  zero. 

The  usual  method  for  specifying  concentrations  of  droplets  at  the  referenc'i 
height  and  the  radii  of  the  corresponding  dry  nuclei  is  as  follows: 

®  User  input  values  for  the  visibility  (meteorological  range)  at  the  reference 
height  and  the  number  of  droplet  size  classes  are  required  in  addition  to  the 
input  specified  above. 

•  The  visibility  and  relative  humidity  are  used  to  define  the  parameters  of  a 
Shettle-Fenn  [3]  bimodal  lognormal  c’istribution  at  the  reference  height. 
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A  lower  dry  particle  (condensation  nucleus)  size  limit  is  set.  The  value  is 
used  to  calculate  the  lower  limit  for  droplet  size  at  the  reference  height. 


•  Because  some  of  the  Shettle-Fenn  distribution  is  not  used,  the  input  value 
of  the  visibility  is  expected  to  be  smaller  than  the  observed  (or  desired) 
meteorological  range  (visibility)  at  thr  reference  height.  To  maintain 
self-consistency,  the  user  is  given  the  option  of  dividing  each  of  the  Ni  and 
N2  parameters  of  the  original  Shettle-Fenn  distribution  by  a  constant.  The 
same  constant  is  typically  chosen  for  use  with  N2  as  with  Nl,  and  the  value 
selected  is  such  that  the  calculated  visibility  (meteorological  range)  at  the 
reference  height  equals  the  input  value  (to  within  a  small  roundoff  erroi). 

•  For  each  size  class,  the  droplet  concentration  and  representative  droplet 
radius  are  calculated  by  integrating  an  appropriate  range  of  the  distribution 
function  and  calculating  the  mean  volume  droplet  radius  for  that  range, 
respectively.  Radii  of  the  dry  nuclei  are  calculated  by  assuming  that  all  of 
the  droplets  at  the  reference  height  are  in  equilibrium  with  water  vapor  in 
air  in  the  cluster  sphere  and  using  an  equilibrium  droplet  radius  equation 
based  on  Hanel  and  Lehmann’s  work  [2]  for  the  type  of  nucleus  selected. 

Sometimes  it  is  desired  to  specify  the  value  for  the  maximum  liquid  water 
content  in  terms  of  a  percent  of  the  quasi-adiabatic  (100  percent)  value  that  is 
obtained  if  all  of  the  other  input  parameters  are  the  same.  This  percent  cannot 
be  specified  directly  as  an  input  value.  However,  the  maximum  liquid  water 
content  depends  strongly  on  the  value  of  the  parameter  alapslin,  which  is  the 
coefficient  of  the  difference  z-zlow  (that  is,  height  minus  reference  height,  both 
expressed  in  units  of  cm)  in  the  cluster  sphere  lapse  rate  equation.  This 
parameter  is  not  requested  interactively;  rather,  it  is  expected  to  be  in  input  file 
zbase_b_c,in.  The  quasi-adiabatic  results  can  be  simulated  by  using  a  value  of 
zero  for  the  alapslin  parameter.  Then  nonzero  values  may  be  used  for  the 
alapslin  parameter  in  subsequent  simulations.  That  parameter  can  be  refined 
iteratively  until  the  maximum  liquid  water  content  is  suficiently  close  to  the 
desired  percent  of  the  corresponding  quasi-adiabatic  value. 
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Examples  of  Input  Files 


In  linis  containing  arrows  in  the  examples,  only  the  input  to  the  left  of  the 
arrows  is  read  by  the  program.  The  arrows  and  the  text  to  the  right  of  the 
arrows  are  included  for  understanding  but  are  ignored  by  the  program.  Some 
of  this  text  refers  to  the  variables  z  the  height  (cm)  above  ground  level,  zlow 
the  reference  height  (cm)  above  ground  level,  and  T  the  te  ^perature  (K)  inside 
the  cluster  sphere. 

NO  fE:  If  standard  input  is  to  be  written  into  a  file  for  use  as  redirected  input, 
it  should  be  noted  that  the  exact  form  of  this  file  depends  on  the  options 
selected.  For  example,  if  a  Shettie-Fenn  distribution  is  chosen  to  define  the 
drop  size  distribution  at  the  initial  (reference)  height,  more  lines  of  input  are 
required  to  provide  the  information  needed  to  evaluate  the  distribution 
parameters  than  if  a  Dutican-Low  distribution  were  chosen.  Also,  if  a 
suggested  value  is  declined,  an  extra  line  is  required  to  provide  the  desired 
input  value. 

The  example  of  standard  input  (collected  into  a  fde  for  use  as  redirected  input) 
includes  two  instances  (pressure  at  reference  height  and  initial  effective  vertical 
velocity)  in  which  the  default  values  are  declined  and  the  entered  values  are 
identical  to  the  defaults.  This  practice  provides  an  explicit  record  of  the  input 
values  used  and  makes  it  easier  for  the  user  to  modify  copies  of  the  standard 
input  file  for  cases  in  which  the  default  values  are  not  the  desired  ones. 

Input  File  Named  zbase_b_c,in  for  Topical  Usage 

zbase,  base  of  INVERSION  in  sphere  lapse  rate  equation  in 
meters 

bhbeqn,  b  .n  extra  term  in  sphere  lapse  rate  equation 
chbeqn,  c  in  extra  term  in  s  phere  lapse  rate  equation 
alapslin,  coeff  of  (z-zlow)  in  sphere  lapse  rate  equation 
clturb,  a  form  factor  to  multiply  KKltur 
cZturb,  a  form  factor  to  multiply  KK2tur 
c3turb,  a  form  factor  to  multiply  KK3tur 
KKltur,  the  turbulent  transfer  coefficient  for  momentum 


5.00E+02  - 


1  9E+04 
2.0E-0/ 
6.5843E-10 
0.0 
0.0 
0.0 
0.0 


0.0  «-  KK2tur,  the  turbulent  transfer  coefficient  for  heat 

0.0  «-  KK3tur,  the  turbulent  transfer  coefficient  for  water  vapor 

mixing  ratio 

2.1.2  Input  File  Named  zbasejbj:,in  with  Priestley-Type  Constants 

S.OOE+02  zbase,  base  of  INVERSION  in  sphere  lapse  rate  equation  in 
meters 

l.OE+04  bhbeqn,  b  in  extra  term  in  sphere  iapse  rate  equation 

2.0E-07  chbeqn,  c  in  extra  term  in  sphere  lapse  rate  equation 

6.5843E-i0  •-  alapslin,  coeff  of  (z-zlow)  in  sphere  lapse  rate  equation 
8  ■*-  clturb,  a  form  factor  to  multiply  KKltur 

8  c2turb,  a  form  factor  to  multiply  KK2tur 

8  ♦-  c3turb,  a  form  factor  to  multiply  KK3tur 

30  ^  KKltur,  the  turbulent  transfer  coefficient  for  momentum 

18  KK2tur,  the  turbulent  transfer  coefficient  for  heat 

30  «-  KKStur,  the  turbulent  transfer  coefficient  for  water  vapor 

mixing  ratio 

2.1.3  Input  File  Named  .haze. defaults 

10000  «=*  plotjnterval  in  units  of  0.1  mm 

100000  **  writ_interval  in  units  of  0.1  mm 

2.1.4  Standard  Input  Collected  Into  a  File  for  Use  as  Redirected  Input 

2  *-  moderate  amount  of  output 

5  ♦-  Shettle-Fenn  distribution  for  droplets  at  reference  height 

2  <-  maritime  air  mass 

^  Blank  line  is  interpreted  as  "yes":  Use  divisors  of  N1 
and  N2. 

0.953773  ♦-  divisor  of  original  Shettle-Fenn  parameter  N1 

0.953773  ♦-  divisor  of  original  Shettle-Fenn  parameter  N2 

4.0R-03  lower  cutoff  (^m)  for  dry  radii 
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4 

n 

0.01 

2 

5.0 

550 

n 

500 

0.1 

n 

20 

a 

0.95 

n 

1000 

n 

288.6827 

n 

288.15 


3 

n 

3 

200 

n 

-1 


drop  concentration  is  inversely  proportional  to  sphere 
volume. 

«-  no  drop  fallout 

entrainment  parameter 
reference  (initial)  height  (m) 
visibility  (km)  at  reference  height 
maximum  height  (m)  for  simulation 
«-  This  is  not  top  of  cloud. 

♦-  height  (m)  of  top  of  cloud 
4-  height  increment  (m) 

Do  not  automatically  accept  default  nucleus  type, 
nucleus  type 

4-  Decline  default  relative  humidity  at  reference  height. 

^  relative  humidity  at  reference  height 
*-  Decline  default  pressure  at  reference  height. 

♦-  pressure  (mbar)  at  reference  height 
*-  Decline  default  reference  height  temperature  in  cluster 
sphere. 

•*-  reference  height  temperature  in  cluster  sphere 
^  Decline  default  reference  height  temperature  of  ambient 
environment. 

reference  height  temperature  of  ambient  environment 
*-  Blank  line  is  interpreted  as  "yes":  Keep  number  of  drop 
sizes. 

^  method  for  determining  effective  vertical  velocity  of 
sphere 

Decline  default  initial  effective  vertical  velocity  of 
sphere. 

■*-  initial  value  (cm/s)  for  effective  vertical  velocity  of 
sphere 

sphere  radius  (cm)  at  reference  height 
dT/dz  in  cluster  sphere  is  not  to  be  treated  as  krown. 
*-  Use  default  droplet  radii  to  define  size  class  boundaries. 
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Blank  line:  Yes,  accept  initial  volume  of  condensed 
matter. 

2.2  Questions  Hiat  Can  Be  Answered  Yes  or  No 

Answer  questions  that  can  be  answered  yes  or  no  by  entering  the  first  letter 
(y  or  n)  of  the  desired  response. 

NOTE:  Avoid  entering  anything,  including  a  blank,  before  the  desired  y  or  n. 
The  first  character  entered,  including  a  blank,  enter,  or  return,  is  considered 
a  yes  answer  unless  the  character  is  n  or  N. 

2.3  In  Case  of  an  Unacceptable  Entry  in  Standard  Input 

A  common  response  to  an  unacceptable  entry  in  standard  input  is 

Unracoonixad  option.  PImm  try  again. 

If  running  the  program  interactively,  enter  revised  input  appropriate  for  the  last 
prompt  that  required  input  and  continue  as  prompted.  However,  if  standard 
input  is  being  supplied  by  a  r^irected  input  file,  the  program  is  not  expected 
to  run  properly  after  the  message  is  encountered.  In  the  latter  case,  output  files 
resulting  from  the  run  in  which  the  error  message  was  generated  should  be 
deleted,  and  the  file  being  used  as  standard  input  should  be  edited  to  change  the 
unacceptable  input  to  acceptable  input.  The  program  should  be  executed  again 
from  the  beginning. 
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3.  Example  Dialogue  Between  MACACASM  and  a  User 


This  section  contains  examples  of  interactive  dialogue  between  MACACASM 
and  a  user.  In  general,  the  examples  are  chronologically  arranged  as  they 
occur  during  an  execution  of  the  program.  As  indicated  in  this  report,  the 
presence  or  absence  of  some  dialogue  depends  on  one  or  more  choices  made 
by  the  user. 


B«tv»Mn  bh«  and  of  this  paragraph  and  tha  and  of  3.29,  taxt  printad  in 
this  f rapraaanta  tuxt  writtan  to  (Standard  output  during  axacution 
of  MAC  ASM.  Unlaaa  atandard  output  ia  radiractad,  tha  text  ia 
axpactad  to  appear  on  tha  terminal  acraan.  Whan  user  choices  result  in 
such  taxt  being  skipped,  as  indicated  by  “Skip  to  ..."  or  similar 
wording,  tha  taxt  does  not  appear  on  tha  terminal  acraan. 

Text  printed  in  this  font  is  for  explanation,  clarification,  etc.,  and  does  not 
appear  on  the  terminal  screen  during  execution  of  the  program. 

3.1  Level  of  Output  Desired 

Indicate  level  of  output  desired: 

(1)  data  files  only, 

(2)  data  filas  plus  some  taxt  (not  much  for  each  drop  size  class),  or 

(3)  data  filas  plus  much  text,  including  values  for  each  size  class. 
Enter  1,  2,  or  3 

The  user’s  choice  greatly  affects  the  amount  of  text  to  be  written  into  output 
file  haze. out,  which  is  an  output  text  file  suitable  for  printing  in  a  132-coIumn 
format.  This  choice  does  not  affect  the  exi?*  'd  content  of  data  files 
suitable  for  use  with  plot  software. 

Text  will  be  written  into  output  file  "haze. out". 

3.2  Method  of  Specifying  Dry  Radii 


Indicate  method  of  specifying  dry  radii; 

(1)  User  specifies  each  dry  radius  and  each  number  of  droplets; 

(2)  User  specifies  each  droplet  radius  at  reference  height  and  each 
number  of  droplets; 
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(4)  Ounean-Lokr  distribution  is  ussd  for  droplot  radii  at  rofsrsno* 
hsigbtf  and  boundarioa  of  droplot  also  intorvals  aro  sot  by  usor  or 
oaloulatodi  or 

(5)  Shottlo-ronn  distribution  is  usod  for  droplot  radii  at  roforonoo 
hoight,  and  boundarioa  of  droplot  sixo  intorvals  aro  sot  by  usor  or 
oaleulatod. 

Tho  rolatirs  huaidity  at  roforonoo  hoight  is  usod  in  calculating  dry 
radii  in  aothods  2,  A,  and  S. 

Intor  tho  dosirod  oothod  (1,  2,  4,  or  5)  •* 

(Method  3  has  been  eliminated.) 

Dry  radii  (equivalent  sphere  radii  of  dry  condensation  nuclei)  are  calculated  or 
specified  at  the  initial  (reference)  height.  The  dry  radii  are  used  during  the 
simulation,  except  when  the  associated  droplets  are  removed  by  gravitational 
settling.  If  method  1  is  selected,  the  dry  radii  can  be  specified  directly  (default 
values  also  available).  If  another  method  is  selected,  a  set  of  droplet  (wet) 
radii  is  calculated  or  specified  at  the  initial  height  and  used  with  the  relative 
humidity  to  determine  the  sec  of  dry  radii  by  turning  that  the  droplets  are 
perfect  spheres  at  their  equilibrium  sizes.  If  method  2  is  selected,  the  initial 
droplet  radii  and  concentradons  can  be  specified  directly  (dehmlt  values  also 
available). 

Methods  4  and  5  define  a  continuous  drop  size  distribution  at  the  initial 
(reference)  height  and  represent  this  distribution  by  a  discrete  drop  size 
distribution  at  that  height.  To  define  the  discrete  distribution,  the  radius  ranges 
of  a  number  of  bins  (intervals)  equal  to  the  desired  number  of  size  classes  are 
specified.  It  is  convenient  to  define  one  lower  cutoff  value  and  let  the  program 
calculate  the  remaining  ranges.  A  representative  radius  is  calculated  for  each 
bin  (range  or  interval).  Because  of  the  importance  of  volume  in  calculating 
liquid  water  content,  the  mean  volume  drop  radius  in  a  given  interval  is  used 
as  the  representative  radiu.<i  for  that  interval.  For  each  size  class,  the  initial 
drop  concentration  and  the  representative  drop  radius  are  determined  by 
integration,  with  respect  to  drop  radius,  using  the  drop  size  distribution 
expressed  as  a  function  of  drop  radius. 

The  continuous  drop  size  distribution  type  used  at  the  initial  (reference)  height 
may  be  different  depending  on  the  method.  If  method  4  is  selected,  a 


Duncan-Low  distribution  is  used  for  droplet  radii  at  the  initial  height.  The 
distribution  is  a  bimodal  gamma  distribution,  parameters  of  which  are  functions 
of  visibility.  If  method  5  is  selected,  a  Shettle-Fenn  distribution  is  used  for 
droplet  radii  at  the  ir  itial  height.  A  Shettie-Fenn  distribution  is  a  bimodal 
lognormal  distribution,  of  which  parameters  are  functions  of  visibility  and 
relative  humidity. 

If  method  1  or  2  is  selected,  skip  to  3.5. 

If  method  4  is  selected,  skip  to  3.4. 

3.3  If  Method  5  is  Selected 

If  method  5  is  selected,  the  following  prompts  appear  in  standard  output  (on 
terminal  screen  if  program  is  being  run  interactively): 

ilpacify  which  Sh«ttl«-Fcnn  mod*!  is  to  ba  usad, 

U)  rural, 

{2)  Mritima,  or 
(3)  urban. 

Bntar  nunbar  of  dasirad  modal  (1,  2,  or  3)  •* 

The  user  is  presented  with  a  choice  of  three  air  mass  types  (maritime,  rural, 
arMl  urban)  and  is  asked  to  select  one.  Near  a  coastline  with  wind  blowing  in 
from  the  ocean,  a  maritime  air  mass  is  the  appropriate  selection.  In  a  city,  an 
urban  air  mass  is  appropriate.  Inland  and  not  near  a  city,  a  rural  air  mass  is 
appropriate. 

Do  you  want  to  divids  at  least  on#  of  ths  Shsttle-Fsnn 
parameters  N1  and  N2  by  a  constant  (yes  or  no) 7 

Answer  y  or  n. 

Sometimes  the  visibility,  calculated  as  (LN  50)/(extinction  coefficient,  km  ',  for 
550  nm  wavelength),  does  not  agree  exactly  with  the  input  value  at  the  initial 
(reference)  height.  Agreement  can  be  forced  by  adjusting  the  original 
Shettle-Fenn  parameters  N1  and  N2,  which  are  used  to  specify  the 
concentrations  of  droplets  associated  with  the  respective  modes.  If  such 
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adjustment  is  necessary,  N1  and  N2  should  be  divided  by  the  same  constant. 
Although  the  capability  exists  for  entering  different  values  for  the  two  divisors, 
unequal  values  do  not  preserve  the  shape  of  the  unaltered  Shettle-Fenn  drop 
size  distribution.  If  the  preceding  question  was  answered  yes,  the  user  is 
prompted  for  values  of  constants  as  follows: 


Entar  positiva  constant  usad  to  divida  N1  •* 
Bntar  positlva  constant  usad  to  divida  N2  -* 

3.4  Lower  Cutoff  for  Dry  Particle  Radius 


Plaasa  dafina  lowar  cutoff t  Sntar  amallast  dry  particla  radius  (in 

micromatars )  to  ba  allowad  (or  0  for  no  cutoff)  -* 

Because  of  the  kelvin  cutoff,  extremely  small  droplets  are  not  expected  to 
become  activated  in  typical  cloud  and  subcloud  conditions.  The  impact  of  such 
small  droplets  on  usual  quantities  of  interest,  such  as  total  liquid  water  content 
and  extinction,  is  expected  to  be  minuscule.  Also,  numerical  difficulties  may 
be  encountered  when  attempting  to  determine  an  equilibrium  radius  (droplet 
radius  when  the  dry  radius  is  known  or  vice  versa)  at  a  given  relative  humidity. 
Other  factors  being  equal,  the  ratio  of  wet  radius  to  dry  radius  at  the  point  of 
critical  supersaturation  tends  to  be  smaller  for  smaller  condensation  nuclei. 
For  extremely  small  nuclei,  the  ratio  could  be  smaller  than  the  ratio  used  in 
this  software  to  define  an  initial  guess  for  the  equilibrium  (droplet  or  dry) 
radius.  With  the  initial  guess  being  on  the  wrong  side  of  the  maximum  of  the 
Kohler  curve,  the  solution  algorithm  might  continue  to  change  trial  radius 
values  in  the  wrong  direction  until  a  trial  radius  value  is  so  huge  or  tiny  that 
the  program  aborts.  Therefore,  a  lower  cutoff  is  used  for  the  dry  particle 
radius.  A  typical  value  used  for  this  cutoff  is  0.004  ^m. 
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3.5  Number  of  Droplets  Definition  Option 


Do  you  want  tha  numbar  of  droplata  par  cubic  santlmatar 
(axcapt  aa  modifiad  by  droplata  falling  out)  to  ba 

(1)  constant, 

(2)  proportional  to  tha  dansity  of  dry  air, 

(3)  raad  from  binary  data  fila  "drop. cone. in" ,  or 
<4)  invarsaly  proportional  to  olustar  aphara  voluma, 

which  dapands  on  danaity  in  clustar  aphara 

(Entar  1,  2,  3,  or  4)? 

The  user’s  choice  here  is  called  the  ‘’Number  of  Droplets  Definition  Option" 
later  in  this  documentation. 

Option  1  is  a  simple  approximation.  It  is  not  strictly  consistent  with  a  fixed 
number  of  droplets  in  an  expanding  sphere. 

Option  2  is  a  reasonable  approximation  of  the  effect  of  expansion  on  the 
droplet  concentration. 

Option  3  gives  the  user  increased  flexibility.  However,  care  should  be  taken 
to  avoid  internal  inconsistencies  when  modeling  with  this  option. 

Option  4  is  appropriate  if  the  droplet  concentration  is  modeled  as  though  the 
droplets  were  evenly  distributed  throughout  the  sphere  but  do  not  pass  through 
the  surface  of  the  sphere.  This  is  the  option  typically  used.  It  is  also  the  only 
option  (of  these  four)  with  which  the  user  also  is  allowed  to  set  the  value  of  an 
entrainment  parameter. 

Because  the  radius  of  the  cluster  sphere  is  modeled  as  being  inversely 
proportional  to  the  cube  root  of  the  partial  density  of  dry  air  in  the  cluster 
sphere,  options  2  and  4  seem  to  represent  very  similar  effects  on  the  number 
of  droplets  per  cubic  centimeter.  Again,  note  that  Option  4  is  also  the  only 
option  (of  these  four)  with  which  the  user  also  is  allowed  to  set  the  value  of  an 
entrainment  parameter. 

Unless  option  3  is  selected,  skip  to  3.7. 
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Option  3  Only 


If  the  user  chooses  Number  of  Droplets  Definition  Option  3,  tlie  file 
drop.conc.in  must  be  available  as  an  input  file  in  the  connected  directory.  The 
file  must  be  unformatted  (binary)  and  sequential,  and  it  must  contain  one 
record  for  every  height  level  used  in  numerical  integration.  (The  first  record 
is  for  the  initial  height,  the  second  record  is  for  the  initial  height  plus  the 
height  increment  entered  by  the  user,  and  so  forth.)  Each  record  must  contain 
the  height  (cm)  above  ground  level  followed  by  the  concentrations 
(droplets/cm^)  of  droplets  in  all  of  the  individual  size  classes  in  order  (starting 
with  size  class  1,  followed  by  size  class  2,  and  so  on).  Each  height  read  from 
this  file  is  compared  with  the  current  height  being  used  in  model  execution. 
If  there  is  disagreement  (allowing  for  reasonable  tolerance),  execution  of  the 
model  program  aborts. 

Drop  Sedimentation 

Do  drops  fall  out  of  cluatar  aphar*  wnan  thair  aattling  valocitiea 
axcaad  tha  ascant  valocity  of  this  sphara? 

Answer  y  or  n. 

Inclusion  of  this  option  was  a  first-order  attempt  to  model  effects  of 
sedimentation,  or  gravitational  settling,  of  drops.  If  the  user  answers  yes,  all 
drops  in  a  size  class  are  modeled  as  falling  out  of  the  cluster  spuere  when  the 
sedimentation  velocity  exceeds  the  ascent  velocity  of  the  sphere.  The  drops  are 
modeled  as  disappearing.  There  is  no  provision  for  modeling  incoming  drops 
that  have  fallen  out  of  preceding  spheres.  This  question  typically  is  answered 
no. 


Unless  the  Number  of  Droplets  Definition  Option  4  was  selected,  skip  to  3.9. 

Entrainment 

The  entrainment  parameter  a  can  be  nonzero  only  if  the  Number  of  Droplets 
Definition  Option  4  has  been  selected. 


BNTRAINMBNT: 

-  If  tha  antrai^  '<iant  paramatar  ALPHA  ia  poaitiva^  it  will  ba  uaad,  and 
tha  Prlaatlay  paramatara  will  ba  «at  to  zaro. 

-  If  ALPHA  ia  zaro,  tha  Priaatlay  paramatara  will  ba  uaad  if  praaant 
in  fila  "zbaaa_b__c,in". 

-  Hagativa  valuaa  ara  not  allowad  for^ ALPHA. 

Bntar  antrainmant  paraoatar  ALPHA. 

Turbulence  effects,  including  bodi  entrainment  and  shedding,  can  be 
represented  by  expressions  like  (and  one  analogous  to)  those  given  by  Priestley. 
[4]  Priestley-type  expressions  or  traditional  entrainment  expressions  involving 
the  entrainment  parameter  a  may  be  used,  but  do  not  use  both.  (The 
entrainment  parameter  a  is  represented  as  a,  in  a  detailed  technical  report  [1] 
describing  this  model.) 

One  or  neither  of  the  two  options  presented  may  be  chosen;  however,  the  two 
options  represent  such  similar,  or  related,  effects  that  they  must  not  be  used 
simultaneously.  A  small  entrainment  effect  typically  is  chosen  by  giving  a  a 
small  positive  value  (such  as  0.01).  Including  a  significant  entrainment  effect 
affects  some  of  the  model  output  and  causes  it  to  be  different  from  what  is 
expected  for  a  well-mixed  system  containing  only  one  regime  throughout.  For 
example,  the  simulated  temperature  lapse  rate  may  be  superadiabatic. 

If  a  positive  value  is  entered  for  a,  this  message  appears: 

Priastlay  p«rain«t«rs  zr*  «at  to  zaro. 

If  zero  is  entered  for  a,  values  of  turbulent  transfer  coefficients  and  form 
factors  for  momentum  and  heat  according  to  Priestley  [4]  and  analogous 
parameters  for  water  vapor  mixing  ratio  can  be  used  if  they  have  already  been 
read  from  file  zbase_b_c,in.  Early  in  the  execution  of  the  program,  there  is 
an  attempt  to  read  from  the  file.  Existence  of  the  file  is  optional;  however,  the 
default  value  of  zero  is  used  for  each  of  the  turbulent  transfer  coefficients  and 
form  factors  if  the  file  does  not  exist.  An  example  of  a  file  containing  nonzero 
values  for  Priestley-type  turbulent  transfer  coefficients  is  given  in  section  2. 1 .2. 
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3.9  Initial  Height 


What  is  tha  initial  haight  in  matara  abova  ground  laval? 

The  user  response  defines  the  height  of  the  center  of  the  cluster  sphere  at  the 
beginning  of  the  simulation.  Typically,  2  m  is  used. 

3.10  Initial  Visibility 

What  ii  tha  visibility  in  kilomatars  at  tha  initial  height? 

If  the  initial  drop  size  distribution  is  defined  by  a  Duncan-Low  or  Shettle-Fenn 
distribution,  the  user-specified  visibility  is  used  in  defining  parameters  of  that 
distribution. 

3.11  Height  Limit  and  Height  of  Top  of  Cloud 

What  is  tha  final  haight  in  raatars  abova  ground  laval? 

It  is  recommended  to  use  a  value  a  little  higher  than  the  desired  top  of  the 
cloud. 

Doas  tha  top  of  tha  cloud  occur  at  this  haight? 

Answer  y  or  n  —  n  if  previous  recommendation  is  followed.  If  the  answer  is 
n,  the  following  prompt  appears. 

Enter  height  of  top  of  cloud  in  metars  abova  ground  level  -* 

Enter  the  height. 

3.12  Height  Increment 


What  is  the  height  increment  in  maters? 

This  is  the  height  step  size  to  be  used  in  model  calculations.  Using  smaller 
step  sizes  requires  more  computation  time  and  leads  to  more  accurate  results. 
The  height  step  size  should  be  such  that  the  intervals  for  writing  output  into 
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text  file  haze. out  (default  interval  —  10  m)  and  plot  data  into  various  files 
(default  interval  =  0.2  m)  are  each  an  exact  multiple  of  the  height  step  size 
(the  height  step  size  multiplied  by  some  positive  integer,  which  could  be  1). 
The  default  values  for  the  intervals  are  overridden  if  file  .haze.defaults  exists 
and  contains  two  integers;  one  integer  in  each  of  two  lines.  The  first  integer 
is  the  interval  (in  units  of  0. 1  mm)  for  writing  plot  data  into  various  files,  and 
the  second  integer  is  the  interval  (in  units  of  0. 1  mm)  for  writing  output  into 
text  file  haze.out.  Considering  the  units  used,  it  is  appropriate  to  divide  both 
of  the  integers  by  10,000  to  compare  with  the  value  entered  via  standard  input 
for  the  height  increment  in  meters.  In  high-accuracy  research  mode,  step  sizes 
of  0.1  m  are  often  used  if  the  top  of  the  cloud  is  set  at  500  m  or  lower,  and 
step  sizes  of  0.2  and  0.25  m  are  often  used  if  the  top  of  the  cloud  is  set  at 
1000  m. 

3.13  Type  of  Nucleus 


Kaap  typa  of  nuclaus  (14)? 

Answer  y  or  n.  If  the  answer  is  y,  the  default  nucleus  indicated  by  the  number 
in  parentheses  is  used.  If  the  answer  is  n,  the  user  is  presented  with  the 
following  list  of  available  choices. 


Entor  type] 

1-5)  not  uaed 

6)  Hanal:  NaCl 

7)  Kanel:  (NH4)2S04 

8)  Hanel  Model  li  rural,  summer  1966,  Mainz 

9)  Hanel:  urban,  Jan  1970,  Mainz 

10)  Hanel  Model  6:  rural 

11)  Mizpeh  Ramon,  Israel;  desert,  normal  conditions 

12)  Mizpeh  Ramon,  Israel;  desert,  sand  storm 

13)  Mainz,  Germany;  industrial 

14)  West  Coast  of  Ireland;  maritime 

15-20)  Deuselbach,  Germany;  background,  combined  with  particles 
from  surrounding  industrial  zones  and  local  sources 


Enter  the  integer  corresponding  to  the  desired  nucleus  type.  The  nucleus  type 
chosen  affects  values  of  parameters  used  in  modeling  microphysics  of  droplet 
growth  and  evaporation.  For  nucleus  types  11  through  20,  subtract  10  to 
obtain  the  corresponding  sample  number  reported  by  Hanel  and  Lehmann.  [2] 
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It  is  suggested  that  the  nucleus  type  selected  be  consistent  with  the  air  mass 
type  selected.  Combinations  of  air  mass  type  and  nucleus  type  that  have  been 
used  include:  nral  air  mass  with  nucleus  type  17,  maritime  air  mass  with 
nucleus  type  14  or  20,  and  urban  air  mass  with  nucleus  type  13. 

3.14  Some  Initial  (Reference)  Height  Values 

R«l«tiva  huiTti'^ity,  total  prasaura,  and  temparatures  of  cluster  sphere 
and  ambler'^  environment  at  rafaranca  height: 

Keep  relative  humidity  (  .98000  ]7 

Answer  y  or  n.  If  the  answer  is  y,  the  number  in  parentheses  is  used.  If  the 
answer  is  n,  the  user  is  prompted  to  enter  the  desired  relative  huuildtty: 

Enter  relative  humidity 

Enter  the  desired  relative  humidity  (using  a  scale  on  which  saturation  is 
exactiy  1)  at  the  initial  (reference)  height. 

Ka<^‘‘p  total  pressure  (  1000.0  millibars)? 

Answer  y  or  n.  If  the  answer  is  y,  the  number  in  parentheses  is  used.  If  the 
answer  is  n,  the  user  is  prompted  to  enter  the  desired  total  pressure  at  the 
initial  (reference)  height: 


Enter  total  pressure  in  millibars  (hectopascals)  -> 

Enter  the  desired  value  for  this  pressure. 

Keep  ciuBi  ".emperature  <  273.16  K  )7 

Answer  y  or  n.  the  answer  is  y,  the  number  in  parentheses  is  used.  If  the 
answer  is  n,  the  Uicr  is  prompted  to  enter  the  desired  cluster  temperature  at  the 
initial  (reference)  height: 

Enter  cluster  temperature  in  degrees  kelvin  -» 

Enter  the  desired  value  for  the  temperatu-s. 


Keep  temperature  of  ambient  environment  (  273.16  K  )? 
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Answer  y  or  n.  If  the  answer  is  y,  the  number  in  parentheses  is  used.  If  the 
answer  is  n,  the  user  is  prompted  to  enter  the  desired  temperature  of  the 
ambient  environment  immediately  outside  of  the  cluster  sphere  at  the  initial 
(reference)  height  by: 

Enter  temperature  of  ambient  environment 
in  degrees  kelvin  -* 

Enter  the  desired  value  for  the  temperature. 

3.15  Number  of  Droplet  Sizes 


Keep  number  of  droplet  sizes  (  80  ) ? 

Answer  y  or  n.  If  the  answer  is  y,  the  number  in  parentheses  is  used.  If  the 
answer  is  n,  the  user  is  prompted  to  enter  the  desired  number  of  droplet  sizes: 

Enter  new  number  -* 

Enter  the  desired  number  of  droplet  sizes.  If  the  input  number  is  too  large  for 
the  allotted  array  sizes,  it  will  be  changed  to  the  maximum  allowed  value,  and 
a  message  will  be  written.  An  example  of  such  a  message  generated  when  the 
maximum  allowed  value  is  80  follows: 


Owing  to  limit  on  number  of  partitions, 
this  quantity  has  been  reset  to  80. 

3.16  Effective  Vertical  Velocity  of  Cluster  Sphere 


choose  method  for  determining  effective  vertical  velocity  of  cluster: 
Effective  vertical  velocity  of  cluster  is  to  be 
(1)  constant  (value  entered  by  user), 

^2)  read  from  file  eff_vert_v.in, 

(3)  calculated  using  cluster  sphere  buoyancy,  or 

(4)  calculated  as  a  linear  function  of  height. 

Enter  1,  2,  3,  or  4 

Method  3  typically  is  used. 

Unless  method  2  is  chosen,  skip  to  3.18. 
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3.17  If  Vertical  Velocity  is  to  be  Read  from  a  File 


If  the  effective  vertical  velocity  is  to  be  read  from  file  eff_vert_v.in 
(method  2),  the  file  must  be  available  in  the  connected  directory.  The  file  is 
expected  to  be  an  ASCII  file  with  each  line  containing  three  floating-point 
numbers.  The  exact  format  used  to  express  each  number  is  not  critical  as  long 
as  the  numbers  are  in  the  expected  sequence  and  can  be  read  properly  using  a 
free-format  read  statement.  The  sequence  for  each  line  is  expected  to  be: 
height  (cm)  above  ground  level,  total  elapsed  time  (s),  and  effective  vertical 
velocity  of  cluster  sphere  (cm/s)  at  the  listed  height.  The  second  number  is 
read  and  ignored;  however,  the  first  and  third  values  are  used,  so  the  numbers 
must  be  in  the  proper  sequence.  The  height  values  must  be  identical  (within 
reasonable  tolerance)  to  the  heights  of  the  levels  used  in  model  calculations 
within  the  program.  If  there  is  a  disagreement  in  the  height  value,  the  program 
aborts  with  the  following  message: 


STOP  Unexpected  height  encountered  in  file  ef f__vert_v. in 

To  ensure  that  the  heights  of  file  eff_vert__v.in  correspond  with  the  numerical 
integration  height  levels,  the  initial  height  given  in  the  file  must  correspond 
with  the  value  entered  interactively  (or  in  a  redirected  standard  input  file)  for 
the  initial  height  in  meters  above  ground  level  (notice  the  units)  when  running 
the  Rachele-Kilmer  microphysics  model  program  (for  example,  200  cm  and 
2  m,  respectively).  Also,  the  height  interval  (difference  in  height  values  in 
consecutive  lines)  in  file  eff_vert_v.in  must  correspond  with  the  value  entered 
interactively  (or  in  a  redirected  standard  input  file)  for  the  height  increment  in 
meters  (for  example,  10  cm  and  0.1  m,  respectively). 

Note  that  the  structure  expected  for  file  eff_vert._v.in  is  exactly  the  same  as  that 
used  for  the  output  file  eff_vert_v:out.  This  relationship  could  help  a 
researcher  study  the  effect  of  varying  only  one  variable.  This  researcher  might 
initially  choose  method  3  for  determining  the  effective  vertical  velocity,  while 
making  sure  that  the  first  number  in  file  .haze. defaults  is  exactly  10,000  limes 
the  value  that  is  to  be  er.tered  interactively  (or  in  a  redirected  standard  input 
file)  for  the  height  increment  in  meters  (for  example,  1000  and  u.l, 
respectively)  in  subsequent  runs.  (The  height  interval  of  outiiut  in  file 
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eff_vert_v:out  and  a  number  of  other  output  data  files  is  determined  by  the  first 
number  in  file  .haze. defaults  if  that  file  exists.  However,  the  height  increment 
for  numerical  integration  is  determined  by  the  user  response  when  prompted 
for  the  height  increment  in  meters.  The  two  intervals  are  not  automatically 
equal.)  The  researcher  could  copy  the  output  file  eff_vert_v:out  to  file 
eff_vert_v.in  in  the  appropriate  directory  and  use  this  copy  with  subsequent 
runs. 

Skip  to  3.19. 

3.18  Initial  or  Constant  Value  for  Effective  Vertical  Velocity 

Keep  effective  vertical  velocity  (  3.0000000  cm/ sec  )7 

Answer  y  or  n.  If  the  answer  is  y,  the  value  in  parentheses  is  used  for  the 
initial  effective  vertical  velocity  of  the  cluster  sphere.  If  the  answer  is  n,  the 
user  is  prompted  to  enter  the  desired  value: 

Enter  initial  value  (in  cm/aec)  for  effective 
vertical  velocity  of  cluster  -* 

for  method  3  or  4  for  determining  effective  vertical  velocity,  or 

Enter  constant  value  (in  cm/aec)  for  effective 
vertical  velocity  of  cluster  -» 

for  method  1 . 

Enter  the  requested  value  for  the  initial  (or  constant)  effective  vertical  velocity 
of  the  cluster  sphere. 

3.19  Initial  Value  for  Radius  of  Cluster  Sphere 


Enter  radius  (in  cm)  of  cluster  sphere  at  reference  height  -* 

This  radius  typically  is  set  equal  to  the  initial  (reference)  height  above  ground 
level.  Note  that  the  units  of  the  input  values  are  different  for  these  two  entries 
(meters  for  initial  height  and  centimeters  for  the  initial  radius  of  cluster 
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sphere).  Typical  values  are  2  m  for  the  initial  height  and  200  cm  for  the 
radius. 

3.20  Vertical  Gradient  of  Temperature  in  Cluster  Sphere 

!■  dT/dz  in  cluster  sphere  to  be  treated  as  known? 

Answer  y  or  n.  N  typically  is  entered,  in  which  case  the  following  message 
is  written: 

Both  dT/dz  and  mass  loading  initially  are  unknown  in  cluster  sphere.  An 
equation  from  thamtodynainics  and  the  ideal  gas  law  is  added  in  order  to 
detsirmine  both  dT/dz  and  d(water  vapor  mixing  ratio) /dz  as  unknowns  in 
this  sphere. 

If  n  is  entered,  skip  to  3.22. 

3.21  If  User  Elects  to  Specify  dT/dz  Profile 

What  method  is  to  be  used  to  determine  dT/dz  in  cluster  sphere, 

(1)  moist  adiabatic  approximation, 

(2)  dT/dz  •  constant  defined  by  user,  or 

(3)  calculation  from  function  dtdzfn? 

Enter  1,  2,  or  3  -» 

3,22.  J  If  Method  1  (Moist  Adiabatic  Approximation)  Is  Selected 


Moist  adiabatic  approximation  will  be  used. 

This  message  confirms  the  selection.  The  value  of  the  vertical  gradient  of 
temperature  in  the  cluster  sphere  is  calculated  automatically  at  each  height 
using  the  moist  adiabatic  approximation. 

This  option  seems  applicable  within  a  cloud  or  layer  of  fog;  however,  it  is  not 
recommended  for  modeling  the  subcloud  region. 

Skip  to  3.22. 
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3.21,2  If  Method  2  (dT/dz  =  Constant  Defined  by  User)  Is  Selected 


Ent«r  value  for  dT/dz  in  cluster  sphere: 

Enter  a  floating-point  number  for  the  constant  value  of  the  vertical  gradient  of 
temperature  in  the  cluster  sphere  (in  units  of  kelvin  per  cm).  A  negative 
number  is  used  to  model  the  usual  situation  of  temperature  decreasing  with 
height. 

The  option  (constant  dT/dz  throughout  entire  height  range)  is  not  recommended 
if  both  subcloud  and  in-cloud  regions  are  to  be  modeled,  as  the  temperature 
lapse  rate  within  the  cloud  is  expected  to  be  di^erent  from  that  beneath  the 
cloud. 

Skip  to  3.22. 

3.21.3  If  Method  3  (Calculation  from  Function  dtdzfn)  Is  Selected 

This  option  allows  the  advanced  user  to  supply  a  customized  subprogram 
defining  the  behavior  of  dT/dz  in  the  cluster  sphere  as  a  function  of  height. 
This  may  be  done  by  replacing  the  Rachele-Kilmer  microphysics  model 
program’s  dtdzfn  function  subprogram  and  recompiling  (with  a  Fortran  77 
compiler)  and  loading  the  entire  program.  The  dtdzfn  function  subprogram 
must  have  the  same  name,  define  a  real*8  function  in  Fortran  77,  and  have 
exactly  one  formal  parameter  that  is  a  real*8  variable  to  represent  height  (cm) 
above  ground  level.  The  values  of  the  formal  parameter  and  the  variables  in 
common  must  not  be  changed  in  the  function  subprogram.  An  example  of  such 
a  function  subprogram  is  given: 


real’s  function  dtdzfn  (zcm) 

*  function  to  calculate  dT/dz  (in  units  of  degrees  kelvin  per  cm) 

*  given  z  (height  above  ground  in  cm) 

*  Thin  could  be  replaced  by  a  user-defined  function  by  replacing 

*  this  function  subprogram  with  an  "include"  statement  used  to 

*  insert  a  file  containing  a  user-defined  function  having  the  same 

*  neune,  type,  and  formal  parameter. 
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implicit  nona 

*  ...  /params/  common  block: 

real’s  ahanel,  bhanel,  denH20,  dryvol,  dtdz, 

&  eta,  etazero,  fsave,  fzero,  kay,  Izero,  maovermO,  rdry, 

&  rhozero,  Rvapor,  sigma,  sigmaw,  spvolHSO,  theta,  tzero 

integer  nuetype 

common  /parama/  ahanel,  bhanel,  denH20,  dryvol,  dtdz, 

&  eta,  etazero,  fsave,  fzero,  kay,  izero,  msovermO,  rdry, 

&  rhozero,  Rvapor,  sigma,  sigmaw,  apvolH20,  theta,  tzero, 

&  nuctyp^j 

*  ...  /paraffl2/  common  block: 

integer  dtdzmeth,  nodmode 

real’s  Cp,  CpliqH20,  Cpvapor,  dlnfdz,  dqdz, 

&  eff_vert__v,  epsilon,  grvcgs,  klwater,  latanthv,  mixrat, 

&  mwK20,  pdry,  pdynes,  pi,  pwater,  rdcqs,  rhoacgs,  spechum, 

&  thatam,  zlow 

common  /param2/  Cp,  CpliqM20,  Cpvapor,  dlnfdz,  dqdz, 

&  aff_vert_v,  epsilon,  grvcgs,  klwater,  latenthv,  mixrat, 

6  mwH20,  pdry,  pdynes,  pi,  pwater,  rdcgs,  rhoacgs,  spechum, 

&  thetam,  zlow,  dtdzmeth,  nodmode 

real’s  zcm 

if  (zcm  .le.  1600.01d+00)  than 
dtdzfn  ■  -0.1082d+00  /  zcm 
else 

dtdzfn  ■  -0.2485d+00  /  zcm 
end  if 

return 

end 

In  this  example,  none  of  the  variables  available  through  labeled  common  are 
used  in  calculating  the  value  of  dtdzfn.  However,  the  values  are  being  made 
available  in  case  the  user  wants  to  use  them  in  the  definition  of  dtdzfn.  A  user 


can  redefine  this  function  by  making  a  copy  of  the  source  code  and  replacing 
the  statements  (in  the  copy)  that  correspond  to  the  following: 


if  (zcm  .!•.  1600.01d‘4'00)  than 
dtdzfn  ■  oQ .  10a2d'f00  /  zcm 

•Isa 

dtdzfn  -  >0.248Sdt00  /  zcm 
and  if 


CAUTION 


Potsntlal  for  Ruining  Rstulta 

Tfia  common  statemonts  and  ths  typo  daclaratlons  and  values 
of  the  variables  in  any  common  statements  must  not  ba 
changadl  Also,  Ihe  type  declaration  (real's)  and  value  ol  the 
fonnal  parameter  zcm  must  not  be  changadl  Making  such 
changes  could  lead  to  mismatches  with  other  program  units. 
Such  mismatches  could  lead  in  tttfn  to  errors  severe  enough  to 
ruin  simulated  data. 


Local  variables  may  be  defined  and  used  if  none  of  their  names  are  being  used 
in  common  statements  or  as  the  formal  parameter.  Any  newly  defined 
variables  should  be  declared  in  type  statements  (such  as  real*8  or  integer). 

See  section  7.2  for  meanings  of  selected  variables  used  in  MACACASM . 

3.22  If  Error  Message  Is  Encountered  at  This  Point 

For  a  certain  type  of  conflict  in  input,  the  program  may  terminate  abruptly  at 
this  point  with  the  following  error  message: 

STOP  OOPSt  Wrong  height  when  trying  to  use  zubroutina  davv_gat 

If  this  error  message  is  not  encountered,  skip  to  3.23. 

If  this  error  message  is  produced,  it  probably  results  from  a  height  value  in  Hie 
eff_vert_v.in  being  different  from  the  current  height  being  used  by  the 
Rachele-Kilmer  microphysics  model  program  after  the  user  has  specified  that 
the  effective  vertical  velocity  is  to  be  read  from  that  file  (method  2). 
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The  heights  in  file  eff_vert_v.in  must  agree  (within  a  small  tolerance)  with  the 
numerical  integration  height  levels.  The  initial  height  given  in  the  file  must 
agree  with  the  value  entered  via  standard  input  for  the  initial  height  in  meters 
above  ground  level  (for  example,  200  cm  in  file  eff_vert_v.in  agrees  with  2  m 
entered  via  interactive  input).  Also,  the  difference  in  height  values  in 
consecutive  lines  in  file  eff_vert_v.in  must  agree  with  the  value  entered  via 
standard  input  for  the  height  increment  in  meters  (for  example,  a  10  cm 
difference  in  the  file  agrees  with  0.1  m  entered  via  standard  input). 

If  file  eff_vert_v.in  is  produced  by  copying  or  renaming  (using  the  mv 
command  in  UNIX)  file  eff_vert_v:out  produced  by  a  previous  run  of  this  same 
program,  then  either  of  two  adjustments  can  be  made  to  assure  agreement. 
The  preferred  adjustment  is  to  supply  the  appropriate  standard  input  value  for 
the  initial  (reference)  height  and  the  appropriate  first  number  in  the 
.haze. defaults  file  while  repeating  the  run  used  to  create  the  first  eff_vert__v:out 
file.  While  repeating  that  run,  choose  method  3  to  determine  the  effective 
vertical  velocity.  Also,  the  value  entered  via  standard  input  for  the  initial 
height  in  meters  above  ground  level  should  be  identical  to  the  corresponding 
value  entered  in  subsequent  runs.  The  first  number  in  file  .haze. defaults 
should  be  exactly  10,000  times  the  value  entered  via  standard  input  for  the 
height  increment  in  meters  (for  example,  10(X)  and  0,1,  respectively)  in 
subsequent  runs.  (The  tolerance  allowed  in  testing  agreement  of  heights  is 
small  -  only  0.5  mm.)  Copy  the  eff_vert_v:out  output  file  to  file  eff_vert_v.in 
in  the  appropriate  directory  and  use  in  subsequent  runs. 

The  other  adjustment  that  can  be  made  is  to  repeat  the  run  that  produced  the 
error  message,  but  with  one  or  two  of  the  input  values  modified  as  necessary. 
Change  one  or  both  of  the  values  entered  via  standard  input  for  the  initial 
height  in  meters  above  ground  level  and  the  height  inclement  in  meters  to 
agree  with  the  data  already  present  in  the  eff_vert_v.in  file. 

If  file  eff_vert_v.in  was  produced  by  a  method  other  than  copying  or  renaming 
(mv  command  in  UNIX)  file  eff_vert_v:out  produced  by  a  previous  run  of  this 
same  program,  determine  that  the  numbers  in  each  line  are  in  the  expected 
order.  The  first  number  for  each  line  is  expected  to  be  height  (cm)  above 
ground  level,  and  the  last  of  the  three  numbers  in  each  line  is  expected  to  be 


the  effective  vertical  velocity  (cm/s)  of  the  cluster  sphere.  (The  second  number 
in  each  line  is  irrelevant,  because  it  is  read  and  ignored.  However,  some 
number  must  be  present  as  a  place  holder  in  the  second  position.  This 
seemingly  strange  requirement  arose  so  that  files  eff_vert_v.in  and 
eff_vert_v:out  would  be  analogous:  a  copy  of  the  latter  file  produced  by  one 
run  of  the  program  may  be  used  with  no  change  as  the  former  file  in  a 
subsequent  run.  The  second  number  in  file  eff_vert_v;out  is  the  simulated  total 
elpased  time  used  by  the  cluster  sphere  to  rise  h'om  its  initial  height  to  its 
current  height.) 

3.23  Informative  Message  Concerning  Expected  Output  Data 

In  typical  use,  it  is  desired  that  the  cluster  sphere  would  stop  rising  before 
reaching  the  user-specified  height.  The  following  message  (with  the 
appropriate  number  in  place  of  549)  is  written  to  inform  the  user  how  many 
simulated  data  points  to  expect  if  the  simulated  cluster  sphere  ascends  the  entire 
distance  to  the  user-specified  final  height: 


If  this  aimul&tion  continues  throughout  the  entire  user-specified 
height  range,  549  plot  data  points  are  to  be  saved  in  each  of  the 
plotit  files  a'.zyy  and  for  each  size  class  in  multiple-si zeclass  files 
hazedsim.seque  and  mode2_ndrops.s. 

Unless  method  1  or  2  was  selected  for  specifying  dry  radii,  skip  to  3.26. 

3.24  If  Method  1  Was  Selected  for  Specifying  Dry  Radii 


The  next  prompt  appears  if  method  1  (user  specifies  each  dry  radius  and  each 
number  of  droplets)  was  selected  for  specifying  dry  radii.  Otherwise,  skip  to 


3.25. 


Enter  each  dry  radius  (in  cm),  or  enter  a  negative 
number  at  the  first  prompt  to  use  default  radiit 

E'or  droplet  size  1,  dry  radius  >■ 

If  a  non-negative  number  is  entered  here,  it  is  used  as  the  dry  radius  for  the 
first  (smallest)  size  class,  and  values  for  all  of  the  other  dry  radii  and  for 
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droplet  concentrations  for  every  one  of  the  size  classes  at  the  initial  height  must 
be  entered. 

If  a  negative  number  is  entered  here,  the  dry  radii  and  initial  droplet 
concentrations  are  defined  using  default  values,  estimated  from  Meppen  SO 
Profile  10  at  a  height  of  approximately  56.98  m  above  ground  level  by 
assuming  that  the  relative  humidity  was  exactly  1.0000  at  that  height.  It  was 
assumed  that  the  temperature  was  271.4402864  K  and  the  pressure  was 
1006.0164  mbar  (hPa)  at  that  height.  If  using  this  set  of  dry  radii,  it  is  best 
to  specify  that  there  are  exactly  17  drop  size  classes,  because  the  set  of  default 
dry  radii  comprise  exactly  17  values. 

Skip  to  3.26. 

3.25  If  Method  2  Was  Selected  for  Specifying  Dry  Radii 

The  next  prompt  appears  if  method  2  (user  specifies  each  droplet  radius  at 
reference  height  and  each  number  of  droplets,  and  the  relative  humidity  at  the 
reference  height  is  to  be  used  in  calculating  dry  radii)  was  selected  for 
specifying  dry  radii.  Otherwise,  skip  to  3.26. 


Entar  oach  droplat  radius  (in  cm)  at  rafaranca  haight,  or  enter 
a  negative  number  at  the  first  prompt  to  usa  default  radii: 

For  droplet  size  1,  droplat  radius  > 

If  a  non-negative  number  is  entered,  it  is  used  as  the  droplet  radius  for  the  first 
(smallest)  size  class  at  the  initial  height,  and  the  user  also  is  expected  to  enter 
values  for  all  of  the  other  droplet  radii  at  the  initial  height  and  for  droplet 
concentrations  for  every  one  of  the  size  classes  at  the  initial  height. 

If  a  negative  number  is  entered  here,  the  initial  droplet  radii  and  initial  droplet 
concentrations  are  defined  using  default  values,  obtained  from  Meppen  80 
Profile  10  at  a  height  of  approximately  2.03675  m  above  ground  level.  If 
using  this  set  of  droplet  radii,  it  is  best  to  specify  that  there  are  exactly  15  drop 
size  classes,  because  this  set  of  default  droplet  radii  comprise  exactly 
15  values. 
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3.26  Some  Informative  Output 


Some  of  the  user-specified  or  accepted  values  are  repeated  in  a  message  such 
as  the  following: 


Hue  typ* 
fo 

Po  (dyn«K/cm**2) 

Tr  (dagrMS  kalvin) 
Tr*  (dagraai  kalvin) 
Haight  incramant (cm) 


20 

0.95000000 

1000000.0 

288.68270 

288.15000 

10.000000 


Unless  method  4  or  5  was  selected  for  specifying  dry  radii,  skip  to  3.28. 


3.27  If  Method  4  or  5  Was  Selected  for  Specifying  Dry  Radii 


Tlie  next  prompt  appears  only  if  method  4  or  5  (Duncan-Low  or  Shettle-Fenn 
distribution  is  used  for  droplet  radii  at  reference  height,  boundaries  of  droplet 
size  intervals  are  set  by  the  user  or  calculated,  and  the  relative  humidity  at  the 
reference  height  is  to  be  used  in  calculating  dry  radii)  was  selected  for 
specifying  dry  radii.  Otherwise,  skip  to  3.28. 

Entar  radii  (in  cm)  that  dafina  boundariaa 

of  tha  siza  clasBas  at  rafaranca  haight,  or  antar  a 

negative  number  at  tha  first  prompt  to  use  default  radii: 

Entar  largest  droplet  radius  (in  cm) 
in  size  class  1  -» 

A  negative  number  typically  is  entered  here.  An  algorithm  in  the  computer 
program  calculates  the  radii  that  define  the  boundaries  of  the  size  classes  at  the 
initial  (reference)  height.  If  a  lower  cutoff  is  used,  the  size  classes,  except  the 
last  (largest  radius)  one,  span  an  equal  range  on  a  log  of  wet  radius  scale.  A 
message,  as  follows,  appears  in  the  standard  output: 


Default  bounds  will  be  used. 

For  size  class  1,  tha  largest  droplet  radius  •  5.591832E-07 
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3.28  Optional  Scaling  of  Initial  Volume  of  Condensed  Matter 


The  user  is  informed  of  the  volume  of  conden5;ed  matter  at  the  initial 
(reference)  height.  Multiplying  the  value  by  one  million  gives  an  acceptable 
estimate  of  the  initial  liquid  water  content  in  grams  per  cubic  meter.  Unless 
the  drop  concentration  has  been  read  from  a  file,  the  user  is  asked  whether  this 
volume  of  condensed  matter  is  acceptable. 

Tha  calculatad  voluma  of  oondanaad  mattar,  uaing  tha  original  numbara 
of  droplata  you  aubmittad,  ia  1.41892075B-09  par  cm‘‘3. 

la  thia  voluma  of  oondanaad  mattar  accaptabla? 

Answer  y  or  n.  If  this  volume  is  not  acceptable,  the  user  is  prompted  to  enter 
the  desired  volume  of  condensed  matter  in  units  of  cubic  centimeter  per  cubic 
centimeter. 


Entar  daairad  voluma  of  oondanaad  mattar  in  om^S  par  cm~3 

Each  original  number  of  droplets  is  multiplied  by  the  constant  (same  for  each 
droplet  size  class)  required  to  change  the  liquid  water  content  to  the  desired 
value.  (Note  that  this  is,  in  effect,  the  same  type  of  adjustment  performed  if 
the  user  specifies  one  common  divisor  of  the  Shettle>Fenn  parameters  N1  and 
N2  while  using  method  S  for  specifying  dry  radii.)  An  example  of  a  message 
to  a  user  who  has  requested  that  the  volume  of  condensed  matter  be 
l.OE-8  cubic  cm  per  cubic  centimeter  at  the  initial  height  follows; 


Each  original  number  of  droplets  has  been  multiplied  by  a  constant 
to  make  the  volume  of  condensed  matter  at  reference  level  « 
l.OOOOOOOOE-08 

3.29  Informative  Output 


The  vertical  gradient  of  the  water  vapor  mixing  ratio  (labeled  dq/dz  in  output) 
is  written  to  standard  output  at  appropriate  intervals.  An  example  of  a  set  of 
the  messages  follows; 

At  z  -  200.,  dq/dz  -  -1.3719224965186E-8 

At  z  »  1200.,  dq/dz  ■  -1.34406377971248-8 

At  z  -  2200.,  dq/dz  «  -1.3183242455061E-8 
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At 

z 

m 

3200., 

dq/dz 

m 

-1.2932076797549E~8 

At 

z 

m 

4200., 

dq/dz 

m 

-1.2691344870049B-8 

At 

z 

m 

5200., 

dq/dz 

m 

-1 . 2467374660478B-8 

At 

z 

m 

6200., 

dq/dz 

m 

-1 . 2271683234876B-8 

At 

z 

m 

7200., 

dq/dz 

m 

-1.2127477499B4E-8 

At 

z 

m 

8200., 

dq/dz 

m 

-1.208742595402B-8 

At 

z 

m 

9200., 

dq/dz 

m 

-1.229S0S1673876B-8 

At 

z 

m 

10200., 

dq/dz 

m 

-1 . 3329802849208B-8 

mm 

Saturation  occur a  at 

or  just  bslow  haight  >■  105.8  m 

At 

z 

m 

11200., 

dq/dz 

m 

-2 . 1727937243866B-8 

At 

z 

m 

12200., 

dq/dz 

m 

-3 . 4338202784332B-8 

At 

z 

m 

13200., 

dq/dz 

m 

-3 . 14431831S7762B-8 

At 

z 

m 

14200., 

dq/dz 

m 

-3.0S696239S8481B-8 

At 

z 

m 

1S200., 

dq/dz 

m 

-3.006S42870718B-8 

At 

z 

m 

16200., 

dq/dz 

m 

-2.965S64127314B-8 

At 

z 

m 

17200., 

dq/dz 

m 

-2.92860237051148-8 

At 

z 

m 

18200., 

dq/dz 

m 

-2.8937751480236B-8 

At 

z 

m 

19200., 

dq/dz 

m 

-2 . 8602475204665B-8 

At 

z 

m 

20200., 

dq/dz 

m 

-2.8275932222981E-8 

At 

z 

m 

21200., 

dq/dz 

m 

-2.7955726437958E-8 

At 

z 

m 

22200., 

dq/dz 

m 

-2 . 7640411384011B-8 

At 

z 

m 

23200., 

dq/dz 

m 

-2 . 7329063865713B-8 

At 

z 

m 

24200., 

dq/dz 

m 

-2.7021067494661E-8 

At 

z 

m 

25200., 

dq/dz 

m 

-2.67159949307SB-8 

At 

z 

m 

26200., 

dq/dz 

m 

-2.6413540388721E-8 

At 

z 

m 

27200., 

dq/dz 

m 

-2.6113478990618S-8 

At 

z 

m 

28200., 

dq/dz 

■ 

-2.581S64146967B-8 

At 

z 

m 

29200., 

dq/dz 

m 

-2.5S1989793S889E-8 

At 

z 

m 

30200., 

dq/dz 

m 

-2.S226147240252B-8 

At 

z 

m 

31200., 

dq/dz 

u 

-2.493430983365E-8 

At 

z 

m 

32200., 

dq/dz 

m 

-2 . 4644323030055E-8 

At 

z 

m 

33200., 

dq/dz 

m 

-2 . 4356137903048E-8 

At 

z 

m 

34200., 

dq/dz 

m 

-2.4069717392609E-8 

At 

z 

m 

35200., 

dq/dz 

m 

-2 . 3785035444601E-8 

At 

z 

m 

36200., 

dq/dz 

m 

-2.3502076995477E-8 

At 

z 

m 

37200., 

dq/dz 

m 

-2.3220839093577B-8 

At 

z 

m 

38200. , 

dq/dz 

m 

-2.2941333256936E-8 

At 

z 

m 

39200. , 

dq/dz 

m 

-2.2663S8971489B-8 

At 

z 

m 

40200., 

dq/dz 

m 

-2 .2387663931864E-8 

At 

z 

m 

41200. , 

dq/dz 

m 

-2 . 2113614908308E-8 

At 

z 

m 

42200. , 

dq/dz 

m 

-2 . 1840607912031E-8 

At 

z 

m 

43200. , 

dq/dz 

m 

-2.15S1052S59254E-8 

At 

z 

m 

44200. , 

dq/dz 

m 

-2 . 1036887118052E-8 

At 

z 

m 

45200., 

dq/dz 

m 

-1.8752436999865E-8 

At 

z 

m 

46200. , 

dq/dz 

m 

-7.7180730178962E-9 

At 

z 

m 

47200. , 

dq/dz 

m 

2.9339963662882B-8 

At 

z 

m 

48200. , 

dq/dz 

m 

1.0714209560853B-7 

At 

z 

m 

49200. , 

dq/dz 

m 

1.9831668970619E-7 

At 

z 

m 

50010. , 

dq/dz 

m 

2.5081328678096E-7 
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END  OF  SlMUIATT'^Nt  sph«r«  i«  no  longor  rifling  at  haight  ■  500.20 

matarfl. 

Changa  in  numbar  of  pointa  in  plot  data: 

Bxpactad:  549  Actual:  499 


4.  Files  and  Standard  Input  and  Output 


During  a  successful  simulation,  a  number  of  output  files  are  written  into  the 
connected  directory.  Many  of  the  output  files  are  written  in  plotit  format.  A 
file  in  plotit  format  is  an  ASCII  file  containing  one  or  more  plot  data  sets.  The 
first  line  of  each  of  the  data  sets  begins  with  three  integers  in  3i5  format.  The 
first  integer  is  the  number  of  data  points  in  the  data  set,  and  the  next  two 
integers  are  the  line  type  and  symbol  type  to  be  used  in  plotting.  The  last 
integer  may  be  followed  by  text  on  the  same  line.  The  first  line  is  followed  by 
one  line  for  each  data  point.  Each  of  the  additional  lines  contains  two  floating 
point  numbers,  which  are  the  values  of  the  X  and  Y  coordinates  for  one  point. 

In  the  following  list,  the  name  of  the  file  is  followed  by  a  b^ief  description  of 
the  conten.  s  of  the  file: 

standard  input  ■  See  2.1.4  for  an  example. 

standard  output 

haze. out  -  an  output  text  file  suitable  for  printing  with  a  132-column  format. 

mode2_tidrops.s,  status  =  NEW,  access  =  SEQUENTIAL, 

form  =  UNFORMATTED  -  binary  file  for  possible  use  by  another  program. 
The  concentration  of  drops  attributed  to  mode  2  are  written  into  this  file  when 
the  initial  drop  size  distribution  is  bimodal.  When  combined,  the  information 
in  this  file  and  file  hazedsim.seque  enables  construction  of  vertical  profiles  of 
drop  size  distributions  for  each  of  the  two  modes  separately,  as  well  as  the 
overall  vertical  drop  size  distribution  profile. 

hazedsim.seque,  status  =  NEW,  access  =  SEQUENTIAL, 

form  =  UNFORMATTED  -  a  binary  output  file  to  which  height,  drop  radii 
and  concentrations,  and  dr/dz  data  (where  r  is  drop  radius  and  z  is  height 
above  ground  level)  are  written.  The  information  in  this  file  is  sufficient  to 
define  a  vertical  drop  size  distribution  profile  and  can  be  used  with  a  program 
such  as  EECSwave  (extinction_etc_coefficients_for_8_wavelengths)  (see  5.2) 
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and  Mie  efficiency  i^iciors  to  calculate  simulated  vertical  profiles  of  extinction, 
backscatter,  scattering,  and  absorption  coefficients. 

hazedsim.bin,  access  =  SEQUENTIAL,  form  =  UNFORJ4ATTED  -  a 
relatively  small  binary  output  file  that  complements  file  hazedsim.seque. 
Included  in  the  values  that  it  contains  are  equivalent  radii  of  the  dry 
condensation  nuclei. 

aarhoa.zyy  -  an  ASCII  output  file,  not  in  plotit  format.  Every  line  of  this  file 
contains  a  pair  of  values  for  height  (cm)  above  ground  level  and  density  of  dry 
air  (in  cgs  units)  inside  the  cluster  sphere. 

eff_vert__v.in,  status  =  OLD  -  file  that  must  be  available  as  an  input  file  if  the 
user  specifies  that  the  effective  vertical  velocity  of  cluster  is  to  be  read  from 
file  eff_vert_v.in.  Each  line  should  contain  three  floating-point  numbers  to  be 
read  as  height  (cm)  above  ground  level,  a  value  to  be  ignored,  and  the  effective 
vertical  velocity  (cm/s)  of  the  cluster. 

eff_vert_v:out  -  an  ASCII  output  file  in  which  each  line  contains  three 
floating-point  values:  height  (cm)  above  ground  level,  elapsed  time  (s)  since 
cluster  started  ascending  from  reference  height,  and  effective  vertical  velocity 
(cm/s)  of  the  cluster. 

rdry-l-ndrops  -  an  ASCII  output  file,  in  which  the  first  line  contains  the  number 
of  size  classes,  a  number  representing  the  method  used  for  calculating  number 
of  drops,  and  the  density  of  dry  air  at  reference  level.  The  rest  of  the  file 
contains  the  dry  radii  and  the  number  of  particles  in  each  size  class  at  the 
reference  level. 

envf.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
relative  humidity  in  the  ambient  environment  versus  height  (m)  above  ground 
level.  The  ambient  environment  is  outside  of  the  cluster  sphere  but  in  its 
immediate  vicinity  and  at  the  same  height.  Ri  lative  humidity  values  written  to 
plot  files  (and  used  within  this  program)  are  on  a  scale  on  which 
saturation  =  exactly  1. 
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aadlnf.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  first  derivative  of  the  natural  logarithm  of  relative  humidity  (in  the  cluster 
sphere)  with  respect  to  height  (cm)  above  ground  level  versus  height  (m)  above 
ground  level. 

aadqdz.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  first  derivadve  of  the  water  vapor  mixing  ratio  (in  the  cluster  sphere)  with 
respect  to  height  (cm)  above  ground  level  versus  height  (m)  above  ground 
level. 

aadtdz.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  first  derivative  of  temperature  (K)  (in  the  cluster  sphere)  with  respect  to 
height  (cm)  above  ground  level  versus  height  (m)  above  ground  level. 

aaf.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for  the 
relative  humidity  (in  the  cluster  sphere)  versus  height  (m)  above  ground  level. 

aaq.zyy  •  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for  the 
mixing  ratio  for  water  vapor  (in  the  cluster  sphere)  versus  height  (m)  above 
ground  level. 

aat.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for  the 
temperature  (K)  (in  the  cluster  sphere)  versus  height  (m)  above  ground  level. 

amwtot.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  total  volume  of  condensed  matter  (cmVcm^,  actually  a  volume  ratio)  (in  the 
cluster  sphere)  versus  height  (m)  above  ground  level.  Although  this  total 
volume  of  condensed  matter  includes  both  water  and  material  from 
condensation  nuclei,  and  the  density  of  droplets  might  not  be  exactly 
1.0000  g  c^^  ^  it  is  often  an  acceptable  approximation  to  multiply  this  total 
volume  by  one  million  and  call  the  result  liquid  water  content  or  mass  loading 
in  units  of  g  m''\  Immediately  following  the  three  integers  in  the  first  line  is 
the  dry  volume,  which  is  the  cm^  that  would  be  occupied  by  the  condensation 
nuclei  from  droplets  in  one  cm^  of  space  in  the  cluster  sphere  at  the  reference 
(initial)  height  if  no  water  were  present. 
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pwater.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  partial  pressure  of  water  vapor  (dyn/cm*)  (in  the  cluster  sphere)  versus 
height  (m)  above  ground  level. 

6t*vT.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  temperature  (K)  in  the  ambient  environment  versus  height  (m)  above 
ground  level. 

.haze.defaults,  status  =  OLD:  See  section  2.1.3  for  an  example.  -  The 
information  in  this  file  (if  it  exists)  is  used  to  set  height  intervals  at  which  plot 
data  and  printed  output  are  to  be  written.  This  file  should  contain  two  integers 
-  one  integer  per  line  -  to  define  the  intervals  at  which  (1)  plot  data  are  written 
into  various  output  files  and  (2)  data  are  written  into  file  haze. out,  which  would 
be  suitable  for  printing.  Both  of  these  values  should  be  in  units  of  0.1  mm. 
If  the  file  does  not  exist,  the  default  values  in  the  program  are  used.  (In  either 
case,  the  values  are  modified  and  used  in  such  a  way  that  both  intervals  are 
multiples  of  the  height  step  size  selected.) 

drop. cone. in,  status  *  OLD,  access  *  SEQUENTIAL, 
form  =  UNFORMATTED  -  needed  as  an  input  file  if  the  user  has  chosen  to 
have  drop  concentrations  read  from  a  binary  data  file  by  choosing  response  3 
(read  from  binary  data  file  drop. cone. in)  to  the  question  "Do  you  want  the 
number  of  droplets  per  cubic  centimeter  (except  as  modified  by  droplets  falling 
out)  to  be  ...".  If  this  input  file  is  needed,  it  must  contain  one  record  for  each 
height  level,  and  the  record  must  contain  the  height  (cm)  above  ground  level 
followed  by  the  number  of  drops  per  cm^  for  each  size  class  in  order  starting 
with  size  class  1.  If  the  height  read  from  this  file  does  not  agree  (within 
reasonable  tolerance)  with  the  height  being  modeled,  the  program  terminates 
abruptly  with  an  error  message. 

envq.zyy  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated  for 
the  specific  humidity  in  the  ambient  environment  versus  height  (m)  above 
ground  level. 


LWC.Knoll  -  an  ASCII  output  file  in  plotit  format  containing  data  simulated 
for  the  volume  of  condensed  matter  (cmVm^  contained  in  drops  having  radii 
between  0.2S  and  23. S  /xm  versus  height  (m)  above  ground  level. 

zbase_b_c,in,  status  =  OLD:  See  sections  2.1.1  and  2.1.2  for  examples.  -  an 
optional  ASCII  input  file.  Its  use  is  strongly  recommended  to  set  the  height 
(m)  above  ground  level  of  the  base  of  the  inversion  layer  equal  to  the  height 
above  ground  level  supplied  for  the  top  of  the  cloud  in  interactive  input.  If  file 
zbase_b_c,in  exists,  it  is  expected  to  contain  4  or  10  values  -  one  value  per 
line.  The  quantities  expected  in  this  file  are,  in  order: 

(1)  zbase:  height  (m)  above  ground  level  of  the  base  of  INVERSION  in 
sphere  lapse  rate  equation 

(2)  bhbeqn:  b  in  extra  term  in  sphere  lapse  rate  equation 

(3)  chbeqn:  c  in  extra  term  in  sphere  lapse  rate  equation 

(4)  alapslin:  coefficient  of  the  difference  z-zlow  in  sphere  lapse  rate 
equation 

-KJhe  file  may  end  here  if  the  default  value  of  zero  is  desired  for  every 
one  of  the  remaining  constants.) 

(5)  clturb:  a  form  factor  to  multiply  KKltur 

(6)  c2turb:  a  form  factor  to  multiply  KK2tur 

(7)  c3turb:  a  form  factor  to  multiply  KK3tur 

(8)  KKltur:  the  turbulent  transfer  coefficient  for  momentum 

(9)  KK2tur:  the  turbulent  transfer  coefficient  for  heat 

(10)  KK3tur:  the  turbulent  transfer  coefficient  for  water  vapor  mixing  ratio 

COMB. file,  access  =*  SEQUENTIAL,  form  =  UNFORMATTED  -  a  binary 
file  intended  for  internal  use  to  store  data  temporarily.  It  remains  available  and 
can  be  used  with  a  copy  of  subroutine  wrfiles,  that  reads  stored  data  from  an 
unformatted  file  COMB. file  and  writes  the  same  data  into  individual  files,  most 
of  which  are  in  plotit  format. 

compare_Tdiffs  -  an  output  text  file  (if  created)  in  the  form  of  a  table 
comparing  three  sets  of  temperature  differences  T  -  T',  where  T  is  the  cluster 
temperature  and  T'  is  the  ambient  air  temperature  at  the  same  height.  The 
table  headings  would  be: 
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Height 

(m) 


Modeled 

T-r 


T  -  T'  Calc  T  -  T'  calc  using 

using  approx  approx  eq  21  with 

eq  21  with  dL/dz  from  only  drops 
full  dL/dz  in  Knollenberg  range 

bounds.dat,  access  =  SEQUENTIAL,  form  =  UNFORMATTED  -  a  relatively 
small  binary  output  file  that  contains  two  records.  These  records  contain: 
Record  1:  nsizes 

Record  2:  (bound(j),  j  =  0,  nsizes) 

where  nsizes  (integer)  is  the  number  of  drop  size  classes,  and  the  bound  array 
(of  real*8)  contains  the  values  used  as  boundaries  of  drop  radii  (cm)  in  the 
partitioning  of  the  initial  drop  size  distribution  into  size  classes. 
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5.  Simulation  of  Vertical  ProHles  of  Extinction,  Backscatter, 
Absorption,  and  Scattering  Coefficients  and  Visibility 


Sometimes  it  is  desired  to  simulate  vertical  profiles  of  extinction,  backscatter, 
absorption,  and  scattering  coefficients  and  visibility.  Multiplying  Mie 
efficiency  factors  for  single  spherical  droplets  of  given  composition  and  size  by 
the  product  of  the  number  concentration  of  droplets  of  that  size  and  the 
geometric  cross-sectional  area  per  droplet  and  summing  over  all  droplet  sizes 
yields  extinction,  backscatter,  absorption,  and  total  scattering  coefficients.  The 
meteorological  range  (visibility)  is  the  natural  logarithm  of  SO,  which  is  about 
3.912,  divided  by  the  extinction  coefficient  for  a  wavelength  of  0.55  /xm. 

Vertical  profiles  of  drop  size  distributions,  simulated  using  MACACASM,  can 
be  used  to  suppi '  the  needed  values  of  number  concentration  of  droplets  and 
geometric  cross-sectional  area  per  droplet.  The  data  can  be  used  with  other 
software  that  uses  or  calculates  the  Mie  efficiency  factors  necessary  to  generate 
the  desired  vertical  profiles.  Such  software  has  been  prepared  for  use  with  an 
HP  9000  computer.  This  software  includes  a  pair  of  programs: 
calc_and_store_Mie_efficiency_factors_for_pure__water  and  EECSwave. 
Calc__and_store_Mie_efficiency_factors_for_pure__water  is  used  to  calculate  and 
store  a  very  large  number  of  Mie  efficiency  factors  for  later  lookup  and  use. 
The  factors  may  be  used  many  times  without  being  recalculated.  EECSwave 
is  used  to  read  the  stored  values  of  Mie  efficiency  factors  and  drop  size 
distribution  data  simulated  using  MACACASM  and  calculate  vertical  profiles 
of  extinction,  backscatter,  absorption,  and  scattering  coefficients  and  visibility. 

5.1  Program  to  Calculate  and  Store  Mie  Efficiency  Factors 

The  source  code  for  calc_and_store_Mie_efficiency_factors_for_pure_water  is 
stored  in  file  nkilmer/sourccs/store_Mie_ef.f  at  ARL  West.  Unless  the  files 
Effic. Factors  and  +Effic.Factors  containing  Mie  effi'  iency  factors  as  functions 
of  droplet  radius  are  available,  use  this  program  to  create  files.  Use  of  this 
program  involves  a  large  amount  of  number  crunching,  and  execution  takes 
such  a  long  time  that  this  program  should  be  run  using  at  or  batch.  Following 
is  an  example  of  batch  when  the  (compiled  and  loaded)  executable  file  has  been 
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stored  as  /home/u22/nkilmer/bin/store_Mie_ef.x.  The  following  example 
shov/s  text  entered  at  HP-UX  command  level;  therefore,  control-D  is  used  to 
indicate  the  end  of  the  batch  job  commands: 

batch 

/home/u22/nkilmer/bin/store_Mie_ef.x 

AD 

Mie  calculations  [S]  are  used  with  the  complex  indices  of  refraction  for  pure 
water  to  generate  a  large  number  of  values  for  efficiency  factors  as  functions 
of  droplet  radius  for  eight  wavelengths.  After  the  efficiency  factors  are 
calculated,  they  are  stored  and  used  when  needed.  This  procedure  avoids 
repeating  time-consuming  Mie  calculations,  which  would  have  to  be  done  if 
changes  in  complex  indices  of  refraction  with  droplet  growth  were  considered. 
The  very  large  increase  in  computational  efficiency  enables  calculating  for 
more  closely  spaced  droplet  sizes  than  would  be  otherwise  feasible. 

Parts  of  the  code  for  the  program  were  copied  from  parts  of  program  AGAUS 
(dated  16-Sept- 1987)  and  converted  to  double  precision.  The  code  used  from 
AGAUS  includes  subroutines  MIEGX  and  ANGLE,  most  of  BLOCKDATA, 
and  selected  fragments  from  the  main  program.  The  resulting  double  precision 
version  of  subroutine  MIEGX  is  used  to  do  the  Mie  calculations. 

Eight  wavelengths  -  0.55,  1.06,  3,  4,  5,  8,  10.6,  and  12  /xm  -  are  used.  The 
real  and  imaginary  parts  of  the  index  of  refraction  of  pure  water  are  used  for 
each  wavelength.  These  values  are  set  automatically  within  the  program.  Mie 
efficiency  factors  are  calculated  for  many  sizes  of  pure  water  droplets.  The 
radii  of  the  droplets  and  the  calculated  Mie  efficiency  factors  for  the  eight 
wavelengths  are  stored  in  two  unformatted  sequential  files  (Effic. Factors  and 
+Effic. Factors)  in  the  connected  directory. 


5.2  Program  to  Calculate  Vertical  Profiles  of  Extinction, 
Scattering,  Backscatter,  and  Absorption  Coefficients  and 
Visibility 

The  source  code  for  EECSwave  has  been  stored  in  file 
~  nkilmer/sources/extin_hp.f  at  ARL  West  and  is  appropriate  for  use  with  HP 
9000  systems.  (The  Cray  version  of  the  source  code  for  this  program  has  been 
stored  as  file  ~nkilmer/sources/extinC8w,seq.f  at  ARL  West.)  The  name 
EECSwave  came  from  "Extinction  Etc  Coefficients  for  8  WAVElengths." 
EECSwave  is  used  to  calculate  and  write  extinction,  scattering,  backscatter,  and 
absorption  coefficients  and  visibility  as  functions  of  height. 

Extinction,  scattering,  backscatter,  and  absorption  efficiency  factors  for  eight 
wavelengths  are  read  as  functions  of  droplet  radius  from  the  Effic.Factors  and 
+Effic.Factors  files.  Drop  size  distribution  data  are  read  from  files  written 
using  the  version  of  MACACASM  that  writes  sequential  binary  files  and  is 
suitable  for  use  on  a  Hewlett-Packard  9000  computer.  The  files  expected  from 
a  previous  execution  of  MACACASM  are  named  aaf.zyy  (an  ASCII  file), 
hazedsim.bin,  hazedsim.seque,  and  bounds.dat  (binary  sequential  files).  These 
files  must  exist  in  the  connected  directory.  (Program  EECSwave  creates  a 
temporary  binary  sequential  file.  The  temporary  file,  cross_sect.bin,  may  be 
deleted  after  successful  execution  of  the  program.) 

If  any  droplets  are  smaller  than  the  sizes  for  which  stored  Mie  efficiency 
factors  exist,  a  double  precision  version  of  subroutine  MIEGX  from  program 
AGAUS  is  used  to  calculate  the  additional  efficiency  factors.  The  subroutine 
is  one  of  the  parts  of  the  code  that  were  copied  from  parts  of  program  AGAUS 
(dated  16-Sept- 1987)  and  converted  to  double  precision  as  was  done  for  the 
program  discussed  in  section  5.1. 

EECSwave  uses  the  same  eight  wavelengths  -  0.55,  1.06,  3,  4,  5,  8,  10.6,  and 
12  ^m  -  for  which  efficiency  factors  were  calculated  and  stored  using  a 
previously  mentioned  program  (discussed  in  section  5.1  and  stored  in  file 
~  nkilmer/sources/store_Mie_ef.f).  The  real  and  imaginary  parts  of  the  index 
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of  refraction  of  pure  water  are  used  for  each  wavelength,  and  the  values  are  set 
automatically  within  both  programs. 

Drop  size  distribution  data  from  MACACASM  are  used  to  calculate  bounds 
and  total  cross-sectional  area  for  each  drop  size  class.  Bounds  and 
cross-sectional  area  data  are  used  with  Mie  efficiency  factors  to  calculate 
extinction,  backscatter,  absorption,  and  scattering  coefficients. 

The  smallest  value  allowed  for  any  drop  radius  is  0.004  /xm.  One  way  to 
assure  that  the  input  data  does  not  include  any  smaller  droplet  radii  is  to 
specify  this  value  as  the  lower  cutoff  value  for  the  dry  particle  radius  v/hen 
using  MACACASM. 

Drop  radii  and  size  class  bounds  at  the  reference  height  are  used  to  calculate 
interpolating  fractions  for  estimating  bounds  for  each  size  class  above  the 
reference  height.  A  simple  scheme  is  used  to  apportion  the  total  geometric 
cross-sectional  area  of  droplets  in  a  given  size  class.  This  scheme  is  analogous 
to  subdividing  a  rectangle  (the  width  of  which  represents  the  range  of  radii  in 
a  size  class,  as  in  a  histogram)  into  skinnier  rectangles.  Calculated  efficiency 
factors  already  exist  for  one  representative  radius  of  each  of  these  skinny 
rectangles.  Having  a  large  number  of  subdivisions  per  size  class  leads  to 
smooth  vertical  profiles  of  extinction,  absorption,  and  backscatter  coefficients. 

Although  cloud  and  subcloud  droplets  contain  small  amounts  of  matter  from 
condensation  nuclei,  the  complex  indices  of  refraction  for  pure  water  are  used 
as  a  very  reasonable  approximation.  Use  of  complex  indices  of  refraction  for 
pure  water  has  a  definite  advantage.  The  same  complex  indices  of  refraction 
can  be  used  for  all  drops  at  all  heights.  Thus,  one  set  of  Mie  calculation 
results  obtained  for  very  closely  spaced  intervals  can  be  stored  and  used 
indefinitely,  enabling  calculation  of  very  closely  spaced  Mie  efficiency  factors 
without  requiring  large  amounts  of  computer  time.  With  the  luxury  of  very 
closely  spaced  Mie  efficiency  factors,  a  simple  interpolation  scheme,  which  is 
analogous  to  subdividing  bars  of  a  histogram  into  many  more  narrow  bars,  is 
adequate.  With  such  closely  spaced  Mie  efficiency  factors,  very  smooth 
simulated  data  are  obtained,  even  when  backscatter  is  involved. 
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5.2.1  Parameter  to  Be  Supplied  to  Main  Program 

The  path  name  of  the  directory  that  is  expected  to  contain  the  Effic. Factors  and 
+Effic.  Factors  files  that  contain  values  for  Mie  efficiency  factors  as  functions 
of  droplet  radii  is  sought  as  the  parameter  supplied  to  the  main  program.  For 
example,  if  the  executable  file  of  this  program  is  file  extin_hp.x  in  the 
connected  directory  and  the  path  name  of  the  directory  containing  values  for 
Mie  efficiency  factors  is  /home/u22/nkilmer/Mie,  an  appropriate  command  to 
invoke  the  program  would  be: 

extin_hp.x  /home/u22/nkilmer/Mie 

If  the  path  name  parameter  is  missing  or  blank  and  the  two  files  are  not  found 
in  the  connected  directory,  the  user  is  prompted  to  enter  the  path  name.  (All 
other  input  and  output  files  are  to  be  in  the  connected  directory.)  If  responding 
to  a  prompt,  do  not  enter  any  leading  blanks  before  the  path  name.  The  path 
name  may  be  the  full  path  name  or  relative  to  the  connected  directory. 

5.2.2  Output  Files 

The  file  names  extinct.zyy,  visib.zyy,  scatter.zyy,  backsca.zyy,  and 
absorpt.zyy  in  the  connected  directory  are  reserved  for  output  files.  The 
program  will  abort  if  any  of  these  files  exist  in  the  connected  directory  when 
the  program  is  started. 

Vertical  profiles  of  simulated  extinction,  scattering,  backscatter,  and  absorption 
coefficients  are  to  be  represented  in  output  files  extinct.zyy,  scatter.zyy, 
backsca.zyy,  and  absorpt.zyy,  respectively.  The  four  files  contain  data  for  all 
eight  of  the  wavelengths  listed.  The  output  file  named  visib.zyy  is  defined  only 
for  a  wavelength  of  0.53  /xm,  because  it  contains  simulated  visibility 
(meteorological  range)  data.  All  output  files  are  in  plotit  format.  Each  data 
set  in  these  files  represents  a  vertical  profile  for  one  wavelength.  Each  line 
that  represents  a  data  point  gives  a  value  of  the  quantity  suggested  by  the  file 
name  (not  including  .zyy)  followed  by  the  height  in  meters.  The  first 
mentioned  value  is  in  units  of  km  for  visibility,  km  ‘  sr  '  for  a  backscatter 
coefficient,  or  km  '  for  an  extinction,  scattering,  or  absorption  coefficient. 
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6.  Other  Selected  Programs  That  Can  Be  Used  With  Output 
from  MACACASM 

The  files  ~nkilmer/sources/hp_maca_dsd.f  and  —  nkilmer/sources/wrt_ 
boundsHP.f  contain  programs  that  can  be  used  to  enable  writing  ASCII  files 
containing  data  that  were  written  to  binary  (unformatted)  files  during  execution 
of  MACACASM.  This  capability  could  be  helpful  if  the  information  is  to  be 
copied  to  a  system  of  a  type  different  from  that  on  which  the  binary  files  were 
created. 

6.1  Rewriting  Drop  Size  Distribution  Data 

Program  get_dsd,  which  is  stored  in  file  -  nkilmer/sources/hp_maca__dsd.f  (for 
an  HP  9000  system),  enables  rewriting  data  from  three  binary  sequential  files 
previously  written  using  MACACASM.  (Similar  programs  stored  as 
files  ~nkilmer/sources/maca_dsd_seq.f  and  nkilmer/sources/maca_dsd 
_10  m.f  could  be  used  on  a  Cray  system.)  One  ASCII  file  (dsd.dat)  is  written 
for  possible  archival  storage  of  simulated  vertical  profiles  of  drop  radius  and 
drop  number  concentration  for  each  size  class.  Drop  radius  versus  z  (height 
above  ground  level)  data  may  be  read  directly  from  this  file  (dsd.dat)  using 
plotit. 

Because  a  very  large  file  could  be  created  if  ASCII  representations  of  all  of  the 
data  were  written  in  their  full  real*8  precision,  the  precision  of  the  numbers 
written  in  ASCII  form  is  reduced  considerably.  Also,  the  units  used  for  drop 
radii  are  changed  from  centimct“Ts  to  micrometers,  and  lines  are  formatted  so 
the  output  file  is  relatively  compact  for  an  ASCII  file. 

Three  unformatted  (binary)  sequential  files  (hazedsim.bin,  hazedsim.seque,  and 
mode2_ndrops.s)  should  have  been  written  using  MACACASM  and  should  be 
in  the  connected  directory.  The  output  ASCII  file  is  to  be  named  dsd.dat,  and 
there  must  not  be  any  other  file  having  that  name  in  the  connected  directory. 
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6.2  Rewriting  Bounds  of  Drop  Size  Classes 


The  program  stored  in  file  ~nkilmcr/sources/wrt_boundsHP.f(foran  HP  9000 
system)  enables  rewriting  data  from  an  unformatted  file  named  bounds.dat, 
which  should  have  been  written  using  MACACASM  and  should  be  in  the 
connected  directory.  The  data  representing  drop  radii  that  define  bounds  of 
drop  size  classes  at  the  reference  (initial)  height  are  written  as  an  ASCII  file 
to  be  named  bounds.ascii.  (A  similar  program  stored  as  file 
~  kilmer/sources/write  _bounds.f  could  be  used  on  a  Cray  system.) 
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7.  Selected  Programming  Notes  for  MACACASM 


This  version  uses  real*8  variables,  which  are  double  precision  on  an  HP  9000 
system.  However,  numerical  sensitivity  in  function  xmeax  is  so  great  that 
some  REAL*  16  precision  is  required.  Function  xmeax  is  now  defined  as  a 
REAL*16  function,  and  there  are  some  REAL*16  variables  in  function  expval 
and  subroutine  iduncan. 

7.1  Files  and  Standard  Input  and  Output 

Integer  variables  are  used  for  unit  numbers  for  files  and  standard  input  and 
output  to  simplify  changing  unit  numbers  if  such  changes  are  needed  for 
portability. 

Some  output  data  are  saved  in  files  for  possible  plotting  or  other  subsequent 
treatment.  Names  of  these  files  are  included  in  this  section. 

The  integer  variable  used  for  the  unit  number  is  followed  by  the  name  of  the 
file  in  the  following  list:  (Section  4  contains  descriptions  of  files.) 

stdin  -  standard  input:  See  2. 1 .4  for  an  example. 

stdout  -  standard  output 

unit26  -  file  haze. out 

unit27  -  file  mode2_ndrops.s,  status  =  NEW,  access  =  SEQUENTIAL,  form 
=  UNFORMATTED 

unit28  -  file  hazedsim.seque,  status  =  NEW,  access  =  SEQUENTIAL, 
form  =  UNFORMATTED 

unit29  -  file  hazedsim.bin,  access  =  SEQUENTIAL, 
form  =  UNFORMATTED 


55 


unit32  -  f!le  aarhoa.zyy 


unit35  -  file  eff_vert_v.in,  status  =  OLD 

unit36  -  file  eff_vert_v:out 

units?  -  file  rdry+ndrops 

unit39  -  file  envf.zyy 

unit41  -  file  aadlnf.zyy 

unit42  -  file  aadqdz.zyy 

unit43  -  file  aadtdz.zyy 

unit44  -  file  aaf.zyy 

unit45  -  file  aaq.zyy 

unit46  -  file  aat.zyy 

unit47  -  file  amwtot.zyy 

unit48  -  file  pwater.zyy 

unit49  -  file  envT.zyy 

unit52  -  file  .haze. defaults,  status  =  OLD:  See  section  2.1. 

unit53  -  file  drop.conc.in,  status  =  OLD,  access  = 
form  =  UNFORMATTED 

unit54  -  file  envq.zyy 

units?  -  file  LWC. Knoll 


3  for  an  example. 
SEQUENTIAL, 
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unit61  -  file  zbase_b_c,in,  status  -  OLD:  See  sections  2.1,1  and  2.1.2  for 
examples. 

unit64- file  COMB.file,  access  =  SEQUENTIAL,  form  =  UNFORMATTED 

unit65  -  file  compareJTdiffs 

unit72  -  file  bounds.dat,  access  =:  SEQUENTIAL, 
form  =  UNFORMATTED 

7.2  Meanings  of  Selected  Variables 

alapslin  -  coefficient  of  the  difference  z-zlow  in  sphere  lapse 

rate  equation 

alphae  -  entrainment  parameter  (called  a  in  dialogue  with 

user) 

-•  If  the  entrainment  parameter  ALPHAE  is  positive,  it 
will  be  used,  and  the  Priestley  parameters  will  be  set 
to  zero. 

If  ALPHAE  is  zero,  the  Priestley  parameters  will  be 
used  if  present  in  file  zbase_b_c,in. 

Negative  values  are  not  allowed  for  ALPHAE. 

aTmTp  a  constant  t  be  used  in  calculation  of  envT,  the 

temperature  in  the  ambient  environment;  This 
constant  and  the  constant  bTmTp  are  defined  and  used 
to  calculate  envT  in  such  a  way  that  the  temperature 
difference  (tzero  -  envT)  is  a  linear  function  of  height 
above  ground  level  and  becomes  zero  at  the  top  of  the 
cloud  (where  z  =  ztop). 

awt  -  the  threshold  water  activity,  which  is  approximated 

(according  to  Hanel  and  Lehmann,  [2])  as  "that 
activity  of  water  at  which  the  volume  of  condensed 
water  is  equal  to  the  volume  of  the  dry  material" 
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awtd 

bhbeqn 

bTmTp 

clturb 

c2turb 

c3turb 

chbeqn 

Cp 

CpUqH20 

Cpvapor 

deltat 

denH20 

direction 


the  threshold  water  activity  divided  by  deltat 

b  in  extra  term,  which  is  effective  near  the  cloud  top, 
in  cluster  sphere  lapse  rate  equation 
a  constant  to  be  used  in  calculation  of  envT,  the 
temperature  in  the  ambient  environment  (See 
aTmTp.) 

a  form  factor  to  multiply  KKltur 

a  form  factor  to  multiply  KK2tur 

a  form  factor  to  multiply  KKStur 

c  in  extra  term,  which  is  effective  near  the  cloud  top, 
in  cluster  sphere  lapse  rate  equation 

heat  capacity  (erg  g  *  K‘‘)  of  dry  air  at  constant 
pressure 

specific  heat  of  liquid  water  (erg  g  *  K'‘) 

heat  capacity  (constant  pressure)  of  water  vapor 
(erg  g  ‘  IC‘) 

the  curvature  correction  when  the  water  activity 
equals  the  threshold  water  activity 

density  (g  cm'^)  of  liquid  water 

rears  variable  that  indicates  the  vertical  direction  of 
motion  the  simulated  cluster  sphere  must  take  to  go 
from  starting  height  to  ending  hei£''t; 

+  1.0  =>  upward 
-1.0  =>  downward 


done 


dqqdz 

dryrad(j) 

dryrmin 

dtdz 


dtdzmeth 


This  variable  should  have  only  the  value  +1.0  when 
set,  because  downward  motion  of  the  cluster  sphere 
is  not  supported  by  this  program. 

-  a  logical  variable  used  to  indicate  whether  simulation 
is  complete.  When  this  variable  equals  .true.,  writing 
to  output  files  will  be  completed,  files  will  be  closed, 
and  execution  will  end. 

-  d(mixing  ratio)/dz  at  the  level  at  which  matrix 
equation  is  solved 

-  dry  radius  (cm)  of  condensation  nucleus  (when 
regarded  as  its  equivalent  sphere)  for  size  class  j 

-  smallest  dry  particle  radius  to  be  allowed  (or  0.0  for 
no  cutoff)  (real's ;  input  by  user  if  not  0.0) 

-  d(temperature)/dz  in  cluster  sphere:  sometimes  at  the 
height  level  at  which  uiatrix  equation  is  solved,  and 
sometimes  the  mean  of  (temperature)/dz  values  at  two 
consecutive  height  levels 

-  method  used  to  determine  dT/dz  (T  is  temperature  in 
cluster  sphere): 

1  -  moist  adiabatic  approximation 

2  -  dT/dz  =  constant  defined  by  user 

3  -  calculation  from  function  dtdzfn 

201  -  both  dT/dz  and  mass  loading  initially  are 

unknown  in  cluster  sphere.  An  equation 
from  thermodynamics  and  the  ideal  gas 
law  is  included  to  determine  dT/dz  and 
d(water  vapor  mixing  ratio)/dz  as 
unknowns  in  this  sphere. 
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dtempdz 


d(temperature)/dz  at  the  height  level  the  matrix 
equation  is  solved.  This  variable  is  never  used  to 
store  the  mean  of  dT/dz  values  at  two  consecutive 
levels,  dtdz  is  used  to  store  the  latter  quantity. 

eff_vert_v  -  effective  vertical  velocity  of  ascending  cluster  (cm/s). 

Sometimes  w  is  used  for  this  quantity. 

envdTdz  -  vertical  gradient  of  temperature  in  ambient 

environment 

envmixrat  -  mixing  ratio  for  water  vapor  in  ambient  environment 

at  same  height  as  cluster  sphere 

envT  -  temperature  (K)  in  ambient  environment  at  same 

height  as  cluster  sphere 

envTr  -  temperature  (K)  in  ambient  environment  at  reference 

(initial)  height 

epsilon  -  ratio  of  the  molecular  mass  of  water  to  the  mean 

molecular  mass  of  dry  air 

eta  -  exponential  mass  increase  coefficient 

etazero  -  value  of  eta  extrapolated  for  an  infinitely  dilute 

aqueous  solution 

fall_out_of_sphere  -  logical  variable  used  to  indicate  whether  drops  fall  out 

of  cluster  sphere  when  their  settling  velocities  exceed 
the  ascent  velocity  of  the  sphere 

fallout(i)  -  logical  variable  used  to  indicate  whether  the  droplets 

in  size  class  i  are  present  (if  .false.)  or  absent  because 
of  gravitational  settling  (if  .true.) 
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first 

iroro 

grvcgs 

KKltur 

KK2tur 

KK3tur 

last_eq_si2e 


latenthv 

lawcoO) 

lawncf 


-  a  logical  variable  used  to  indicate  whether  the 
simulated  cluster  sphere  is  still  at  its  first  (reference) 
height  (if  .true.)  or  at  some  other  height  (if  .false.) 

-  relative  humidity  (1.0  -  saturation)  in  cluster  sphere 

-  acceleration  due  to  gravity  (cm  sec'^) 

-  the  turbulent  transfer  coefficient  for  momentum 

-  the  turbulent  transfer  coefficient  for  heat 

-  the  turbulent  transfer  coefficient  for  water  vapor 
mixing  ratio 

-  last  size  class  where  droplets  are  always  at  their 
equilibrium  sizes.  The  variable  can  be  used  if  drops 
in  a  certain  number  of  the  smaller  size  classes  are 
modeled  as  always  being  at  their  equilibrium  sizes. 
In  normal  use  of  the  program,  all  droplets  exhibit 
diffiision-controlled  growth  and  evaporation  above  the 
reference  height;  therefore,  last_eq_size  is  set  equal 
to  0. 

-  specific  enthalpy  of  vaporization  of  water  (erg  g  ‘) 

-  the  coefficient  of  it*  (v\..cre  rr  is  the  ratio  of  dry 
radius  to  wet  radius)  in  the  fitting  polynomial  used  to 
calculate  eta/etazero  as  part  of  the  calculation  of  the 
activity  of  water  in  a  droplet 

-  degree  of  fitting  polynomial  (in  powers  of  the  ratio  of 
dry  radius  to  wet  radius)  used  to  calculate  eta/etazero 
as  part  of  the  calculation  of  the  activity  of  water  in  a 
droplet 
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level 


number  of  current  level  (reference  level  =  1,  etc.) 


meth_for__evv 


mixrat 
mwH20 
Ni  divisor 

N2divisor 

ndaiTay(j) 

nodmode 


store  -  method  for  determining  effective  vertical  velocity  of 
cluster  sphere  (integer  chosen  by  user): 

Effective  vertical  velocity  of  cluster  is 

1  -  constant  (value  entered  by  user) 

2  -  read  from  Hie  eff_vert_v.in 

3  -  calculated  using  cluster  sphere  buoyancy 

4  -  calculated  as  a  linear  function  of  height 

-  mixing  ratio  for  water  vapor  in  cluster  sphere 

-  molecular  weight  of  water 

-  constant  used  to  divide  Shettle-Fenn  parameter  Nl 
(real*8;  input  by  user  if  different  from  1 .0) 

•  constant  used  to  divide  Shettle-Fenn  parameter  N2 
(real*8;  input  by  user  if  different  from  1.0) 

-  number  of  droplets  per  cm^  in  size  class  j 

-  Number  Of  Drops  MODE:  integer  input  by  user 
(integer  used  internally).  The  number  of  droplets/cm^ 
(except  possibly  as  modified  by  droplets  falling  out) 
are 

1(1)  -  constant 

2(2)  -  proportional  to  the  density  of  dry  air 
3(6)  -  read  from  binary  data  file  drop.conc.in 
4(8)  -  inversely  proportional  to  cluster  sphere 
volume,  which  depends  on  density  in 
cluster  sphere 


nsizes 


number  of  droplet  size  classes 


nuctype 

outplevel 


pdry 

pdryref 

pdynes 

plotjnterval 


Pradius 

Pradref 

prelim 

pwater 


nucleus  type  (integer  input  by  user) 
level  of  output  desired: 

1  -  data  files  only 

2  -  data  files  plus  some  text  (not  much  for  each  drop 

size  class) 

3  -  data  files  plus  much  text,  including  values  for 

each  size  class 

partial  pressure  (dyn  cm'^)  of  dry  air 

partial  pressure  (dyn  cm*^)  of  dry  air  at  reference 
(initial)  height 

total  pressure  (dyn  cm*^) 

height  interval  at  which  plot  data  are  written 
(Although  initially  in  units  of  0.1  mm,  integer  is 
converted  to  the  exact  number  of  height  steps  in  the 
interval.) 

radius  (cm)  of  cluster  sphere 

radius  (cm)  of  cluster  sphere  at  reference  height 
a  logical  variable  used  to  indicate  whether  execution 
is  in  preliminary  mode  (.true,  if  in  preliminary  mode; 
.false,  otherwise) 

partial  pressure  (dyn  cm’*)  of  water  vapor 
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pzero 

radius(i) 

rdcgs 

regimcO) 

rhoacgs 

rhoaref 

Rvapor 

SFmodel 


sigmaw 

sizmode 


-  total  pressure  (mbar)  (Although  the  total  pressure  at 
the  initial  height  is  accepted  from  the  user  in  these 
units,  the  corresponding  cgs  variable,  pdynes,  is 
defined  and  used  internally.) 

-  radius  (cm)  of  drop  in  size  class  i 

-  specific  gas  constant  (erg  g'‘  K'*)  for  dry  air 

-  regime  of  curvature  correction  for  size  class  j  (In 
some  cases,  different  equations  are  appropriate  for 
different  regimes.  The  regime  array  is  used  to 
communicate  which  curvature  correction  regime  is  in 
effect  for  each  droplet  size  class.) 

-  density  of  dry  air  (g/cm’)  in  cluster  sphere 

-  density  of  dry  air  (g/cm^)  in  cluster  sphere  at 
reference  (initial)  height 

-  specific  gas  constant  for  water  vapor  (erg  g“‘  K"‘) 

-  an  integer  that  indicates  which  Shettle-Fenn  model  is 
to  be  used  (input  by  user): 

1  -  rural 

2  -  maritime 

3  -  urban 

-  surface  tension  of  liquid  water  at  273. 16  K 

-  method  of  specifying  dry  radii: 


64 


spechum 

spvolH20 

thetam 

tzero 

visib 

wplot 

writ 

writ  interval 


z 


1  -  User  specifies  each  dry  radius  and  each  nunnber 

of  droplets 

2  -  User  specifies  each  droplet  radius  at  reference 

(initial)  height  and  each  number  of  droplets 

4  -  Duncan-Low  distribution  used  for  droplet  radii  at 

reference  height,  and  boundaries  of  droplet  size 
intervals  set  by  user  or  calculated 

5  -  Shettle-Fenn  distribution  used  for  droplet  radii  at 

reference  height,  and  boundaries  of  droplet  size 
intervals  set  by  user  or  calculated 
The  relative  humidity  at  reference  height  is  used  in 
calculating  dry  radii  in  methods  2,  4,  and  S. 

-  specific  humidity 

-  specific  volume  (cm^)  of  pure  liquid  water 

-  damp  potential  temperature 

-  temperature  (K)  in  cluster  sphere 

-  visibility  (km)  at  initial  (reference)  height 

-  a  logical  variable  used  to  indicate  whether  current 
data  are  to  be  written  as  plot  data 

-  a  logical  variable  used  to  indicate  whether  current 
data  are  to  be  written  as  text,  usually  into  file 
haze.out 

-  height  interval  at  which  printed  output  is  written  as 
text,  usually  into  file  haze.out  (Although  initially  in 
units  of  0.1  mm,  integer  is  converted  to  be  the  exact 
number  of  height  steps  in  interval.) 

-  height  (cm)  above  ground  level 
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zbase  •  height  (cm)  above  ground  level  of  the  base  of 

inversion  in  sphere  lapse  rate  equation:  It  is 
recommended  that  the  height  be  identical  to  ztop. 
(Note:  If  file  zbase__b_c,in  exists,  then  zbase  is  read 
from  that  file  in  units  of  meters  and  converted 
immediately  to  units  of  cm.) 

zhigh  -  user-specified  final  height  (cm)  above  ground  level  for 

the  simulated  cluster  sphere.  (Note:  If  simulated 
cluster  sphere  stops  rising  before  this  height  is 
attained,  the  last  height  level  at  which  this  cluster 
sphere  is  rising  would  be  used  as  the  final  height 
modeled.)  (Also  note:  User  input  values  for  zhigh, 
zinc,  zlow,  and  ztop  are  read  in  units  of  meters  and 
converted  to  units  of  cm.) 

zinc  -  height  increment  (step  size)  (cm) 

zlow  -  initial  height  (cm)  above  ground  level  for  the 

simulated  cluster  sphere 

ztop  -  height  of  top  of  cloud  (cm)  above  ground  level:  It  is 

recommended  that  this  neight  be  identical  to  zbase. 

Threshold  values  related  to  curvature  correction  include  awt,  awtd,  and  deltat, 
which  are  in  common  block  /param3/.  Definitions  of  these  variables  are 
included  in  the  preceding  list. 

7.3  General  Notes 

7.3.1  Frequency  of  Output 

If  file  .haze. defaults  exists  in  the  connected  directory,  the  information  in  that 
file  is  used  to  set  height  intervals  at  which  plot  data  and  printed  output  are  to 
be  written.  If  the  file  is  not  present,  default  values  are  used  to  set  the 
intervals.  The  default  height  interval  at  which  plot  data  are  written  to  various 
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files  is  0.2  m.  The  default  initial  value  for  the  associated  variable, 
plotjnterval,  is  2000  (units  of  0.1  mm).  The  default  height  interval  at  which 
printed  output  is  written  to  file  ha^.out,  which  would  be  suitable  for  printing, 
is  10  m.  The  default  initial  value  for  the  associated  variable,  writ_interval,  is 
100000  (units  of  0. 1  mm). 

The  two  associated  variables  are  modified  in  such  a  way  that  the  plot  and  write 
intervals  will  be  exact  multiples  of  the  height  step  size  selected.  The  modified 
definition  of  each  of  these  associated  variables  (plot_interval  and  writ_interval) 
will  be  the  number  of  height  steps  (of  length  zinc)  in  the  interval. 

7.3.2  Default  Values  Used  if  File  "7jbase_b_c,in“  Is  Not  Present 

If  file  zbase_b_c,in  is  not  present  in  the  connected  directory,  default  values  are 
used  in  place  of  the  constants  that  would  have  been  read  from  the  file.  The 
default  values  include  these  values  for  constants  in  the  cluster  sphere  lapse  rate 
equation: 

zbase  -  2.55d+04  for  base  of  inversion  (cm) 

bhbeqn  =  O.d+00 

chbeqn  =  4.d-04 

alapslin  =  O.d+00 

The  default  values  effectively  eliminate  the  terms  in  v'hich  bhbeqn,  chbeqn, 
and  alapslin  occur. 

Default  values  also  include  these  values  for  Priestley-type  form  factors  and 
turbulent  transfer  coefficients  at  the  reference  height: 

clturb  =  O.d-l-00 

c2turb  =  O.d-l-00 

c3turb  =  0.d+(X) 

KKltur  =  O.d-hOO 

KK2tur  =  O.d+00 

KK3tur  =  O.d+OO 
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7.3.3  Specific  Enthalpy  of  Vaporization  of  Water 

Specitic  enthalpy  of  vaporization  of  water  (in  units  of  erg  g'O  is  based  on  data 
given  by  Haar  et  al.  [6]  and  is  calculated  according  to  this  equation: 

latenthv  =  2.5005d+ 10  -  2.36d+07  *  (tzero  -  273.16d+00) 

7.3.4  Treatment  of  Mass  Due  to  Nuclei 

The  density  of  liquid  water  (denH20)  also  i  used  as  the  density  of  droplets  in 
calculating  the  liquid  water  concentration,  although  the  droplets  are  not 
composed  of  pure  water.  The  mass  of  the  nuclei  is  included  in  calculating  the 
liquid  water  content  but  is  not  included  in  calculating  Cpstar  and  the  liquid 
water  mixing  ratio  (a  variable  named  water)  in  the  cluster  sphere. 

7.3.5  Determination  of  Size  Distributions 

The  simulated  dry  particle  size  (radius)  distribution  and/or  the  simulated  drop 
size  (wet  radius)  distribution  inside  a  representative  cluster  sphere  are 
determined  at  the  reference  (initial)  height.  Either  (sizmode  1)  the  dry  particle 
size  (radius)  distribution  is  to  be  specified  directly,  or  (sizmode  2,  4,  or  S)  the 
drop  size  (wet  radius)  distribution  is  to  be  specified  and  the  dry  particle  size 
distribution  is  to  be  calculated  assuming  (for  the  reference  height  only)  that  the 
system  s  at  equilibrium  inside  the  cluster  sphere. 

If  sizmode  is  4  or  5,  parameters  are  determined  for  the  appropriate  continuous 
drop  size  distribution,  which  is  later  partitioned  into  a  discrete  distribution  with 
the  requested  (or  default)  number  of  size  classes.  After  being  determined, 
discrete  drop  sizes  will  be  used  with  other  reference  height  information 
(including  relative  humidity)  and  constants  associated  with  the  selected 
condensation  nucleus  type  to  calculate  the  dry  particle  size  distribution. 

When  a  Shettle-Fenn  or  Duncan-Low  drop  size  distribution  is  selected,  the  high 
end  of  that  distribution  is  retained.  The  lower  cutoff  corresponds  to  a  dry 
particle  radius  specified  by  the  user. 
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7.3.6  Ensuring  Valid  Values  at  Reference  Height 

To  ensure  calculating  valid  values  at  the  reference  height,  the  first  complete 
execution  of  the  predictor-corrector  loop  is  performed  without  changing  height. 
Flags  are  set  to  indicate  that  execution  is  in  preliminary  mode  and  to  suppress 
storing  plot  data  and  printing  while  in  preliminary  mode. 

7.3.7  Drop  Sedimentation 

If  the  user  has  answered  yes  to  the  question  "Do  drops  fall  out  of  cluster 
sphere  when  their  settling  velocities  exceed  the  ascent  velocity  of  this  sphere?", 
the  sedimentation  velocities  of  drops  are  compared  against  the  ascent  velocity 
of  the  cluster  sphere.  Any  droplet  size  class  in  which  drops  have  sedimentation 
velocities  greater  than  the  ascent  velocity  of  the  cluster  sphere  are  eliminated 
and  are  not  restored  for  the  remainder  of  the  simulation.  If  any  eliminated  size 
class  has  a  smaller  index  (sizeclass  or  i)  than  any  retained  size  class,  size 
classes  still  present  with  higher  indices  also  are  eliminated  and  a  warning 
message  is  written.  The  program  is  terminated  if  all  size  classes  are 
eliminated. 

7.3.8  Assuring  That  Final  Height  Conforms  to  User*s  Request 

The  last  height  is  to  be  exactly  what  the  user  has  specified,  although  it  might 
require  using  a  smaller  increment  (zinc)  for  the  final  step. 

7.3.9  Temporary  Data  Storage 

Data  are  stored  temporarily  in  a  combined  unformatted  file  so  a  number  of 
individual  ASCII  files  can  be  written  in  plotit  format,  which  requires  having 
the  number  of  data  points  available  when  the  first  line  is  written.  Near  the  end 
of  execution  of  the  program,  data  stored  in  the  combined  unformatted  file  are 
read  and  written  into  individual  ASCII  files. 
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7.3,10  Temperature  in  the  Ambient  Environment 


Calculation  of  temperature  in  the  ambient  environment  is  based  on  specifying 
that  the  difference  between  temperatures  in  the  cluster  sphere  and  the  ambient 
environment  varies  linearly  with  height  from  the  value  at  the  initial  (reference) 
height  to  zero  at  the  top  of  the  cloud. 

7.3.11  Turbulence  Effects 

Turbulence  effects  may  be  represented  either  by  traditional  entrainment 
expressions  or  by  expressions  like  (and  one  analogous  to)  expressions  given  by 
Priestley.  [4]  Use  the  Priestley-type  expressions  or  the  traditional  entrainment 
expressions  involving  the  constant  alphae,  but  do  not  use  both. 

7.3.12  Vertical  Acceleration  of  Cluster  Sphere 

The  equation  for  vertical  acceleration  of  the  cluster  sphere  does  not  include  a 
Coriolis  force  term,  that  would  depend  on  the  east-west  component  of  wind 
velocity.  One  possible  method  of  estimating  the  east- west  velocity  component 
uses  similarity  expressions  above  the  reference  height.  However,  it  is  intended 
that  similarity  theory  not  be  used  at  all  in  this  model. 

7.3.13  Radii  of  Droplets  and  Particles 

Radii  are  calculated  as  though  droplets  and  dry  nuclei  were  perfect  spheres. 
The  dryrad  array  stores  radii  of  dry  nuclei.  The  simple  variable  rdry  is  used 
for  the  dry  nucleus  of  the  current  size  class.  The  radius  array  stores  radii  of 
droplets.  The  simple  variable  rwet  sometimes  represents  the  radius  of  a  droplet 
in  the  current  size  class  and  a  related  quantity  (such  as  the  droplet  radius  if  it 
were  in  equilibrium). 

The  initial  guess  used  when  solving  iteratively  for  a  droplet  radius  is  1 .25  times 
the  dry  radius  in  subroutine  solver,  which  is  used  to  calculate  equilibrium 
droplet  radii  at  the  reference  height.  The  initial  guess  is  not  always 
appropriate.  The  program  may  terminate  prematurely  if  the  lower  cutoff  dry 
radius  is  set  to  zero  or  a  very  small  value. 
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A  nonequilibrium  dr/dz  equation  from  Henry  Rachele,  [7]  revised  to 
incorporate  modified  hydrostatic  approximation,  is  used  for  all  droplet  size 
classes  for  all  heights  above  the  reference  height.  There  is  no  automatic 
changing  of  any  size  class  from  equilibrium  to  nonequilibrium,  except  that  all 
size  classes  are  considered  to  be  in  equilibrium  at  the  reference  height,  and  all 
size  classes  are  considered  to  obey  the  nonequilibrium  dr/dz  equation  above  the 
reference  height. 

7.3A4  Units 

Centimeter-gram-second  (cgs)  units  are  used  in  this  program.  Pressures 
(represented  by  variables  named  pdry,  pdynes,  pwater,  etc.)  used  in  typical 
calculations  are  in  units  of  dyn  cm'^. 

7.3  J5  Hydrostatic  Approximation 

Hydrostatic  approximation  is  based  on  values  in  the  ambient  environment 
outside  of  the  cluster  sphere  but  at  the  same  height  as  the  cluster  sphere.  No 
condensed  matter  is  considered,  and  the  variable  envmixrat  (instead  of  mixrat) 
influences  the  hydrostatic  approximation. 

7.4  Selected  Subprograms 

7.4.1  Function  curvcor 

Function  curvcor  calculates  curvature  correction  and  some  related  values. 

7.4.2  Subroutine  dcurvco 

Subroutine  dcurvco  calculates  partial  derivatives  of  the  natural  logarithm  of 
curvature  correction.  Subroutine  rhiunc  should  be  called  before  subroutine 
dcurvco  is  called,  because  the  value  of  awtd  (the  threshold  water  activity 
divided  by  deltat)  that  might  be  needed  in  subroutine  dcurvco,  is  calculated  in 
curvcor,  a  function  that  is  used  in  rhfiinc. 
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7,4,3 


Subroutine  dew _put 


Subroutine  dew_put  stores  reference  height  values  for  effective  vertical 
velocity  and  its  derivative  with  respect  to  height.  If  metli_for_cvv_storc  is  2, 
derivative  of  effective  vertical  velocity  is  stored  for  each  height,  although  the 
values  for  the  reference  height  and  final  height  actually  are  calculated  for  the 
reference  height  plus  one-half  of  a  height  increment  and  the  final  height  minus 
one-half  of  a  height  increment,  respectively. 

7.4.4  Subroutine  DROGESUP 

Subroutine  DROGESUP  (DROplet  Growth  Equation  SetUP)  calculates  the 
coefficients  (dqdz_coeff  and  dtdz_coeff)  of  the  vertical  gradients  of  mixing 
ratio  and  temperature  and  the  right  hand  side  (rhs)  of  an  equation  involving  the 
vertical  gradient  of  the  radius  of  a  droplet  of  size  class  sizeclass.  These  values 
may  be  used  in  defining  a  row,  the  position  of  which  is  equal  to  the  value  of 
sizeclass,  of  a  matrix  equation  when  droplets  in  this  size  class  exhibit 
diffusion-controlled  (nonequilibrium)  growth. 

7.4.5  Subroutine  environ 

Subroutine  environ  calculates  values  in  the  ambient  environment.  It  is  assumed 
that  (1)  no  condensed  matter  is  in  the  ambient  environment  and  (2)  the  total 
pressure  (dry  air  plus  water  vapor)  in  the  ambient  environment  is  identical  to 
the  total  pressure  in  the  ascending  cluster  sphere.  The  relative  humidity  in  the 
ambient  environment  is  set  equal  to  the  relative  humidity  in  that  cluster  sphere. 

7.4.6  Subroutine  errfun 

Subroutine  errfun  calculates  the  error  function  and  the  complement  of  the  error 
function  according  to  Cody.  [8]  Its  formal  parameters  are: 

X  -  argument  for  erf  or  erfc  (real*8;  input) 
erf  -  value  of  error  function  (rear8;  output) 

erfc  -  value  of  complement  of  error  function  (real*8;  output) 
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which  -  indicates  which  result  is  calculated  directly  from  a  rational  function 
(integer;  output): 

1  if  erf  is  calculated  directly  (and  erfc  from  erf) 

2  if  erfc  is  calculated  directly  (and  erf  from  erfc) 

7.4.7  Subroutine  evvjcalc 

Subroutine  ew_calc  calculates  the  sphere  radius  (Pradius)  using  the  ratio  of 
partial  density  of  dry  air  (in  the  sphere)  to  that  at  the  reference  height.  This 
calculation  does  not  count  any  of  the  entrained  matter  as  part  of  the  sphere. 

If  appropriate  (that  is,  if  meth_for_ew_store  is  3  or  4),  the  sphere  radius  and 
other  quantities  are  used  to  calculate  dewdz  (i.e.,  dw/dz),  which  is  used  to 
calculate  w  by  numerical  integration. 

The  formal  parameters  for  subroutine  ew_calc  are: 

z  -  height  (cm)  (input) 

w  -  effective  vertical  velocity  of  cluster  (without  including  any  turbulent 
fluctuations)  (ouq)ut;  value  is  also  put  into  labeled  common) 

7.4.8  Function  expval 

The  function  expval  (with  formal  parameters  n,  rlcm,  and  r2cm)  determines 
the  expectation  value  of  r*  for  the  interval  of  r  from  rlcm  to  r2cm.  The 
method  includes  evaluating  closed-form  solutions  of  sums  of  definite  integrals 
of  the  Duncan-Low  n(r)  density  function. 

The  formal  parameters  rlcm  and  r2cm  are  droplet  radii  in  units  of  cm. 
Although  these  are  passed  unchanged  to  subroutine  iduncan,  the  equivalent 
radii  (rl  and  r2)  in  micrometers  are  used  when  passing  values  to  function 
xmeax. 
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7,4,9  Subroutine  getdrs 


Subroutine  getdrs  sets  up  a  matrix  equation  that  is  solved  for  d'*/dz  for  each 
size  class,  (1/f)  df/dz  sometimes  (if  at  reference  height  or  if  entrainment  is  not 
considered),  and  dWv/dz  and/or  dT/dz  in  cluster  sphere.  (Here,  f  represents 
relative  humidity  on  a  scale  on  which  saturation  is  exactly  1,  Wv  represents 
water  vapor  mixing  rado  in  the  cluster  sphere,  and  T  represents  temperature 
(IQ  in  the  cluster  sphere.)  Subroutine  PADLIMES  does  the  solving. 
Calculated  values  are  stored  in  the  drdz  array  and  the  simple  variables  dlnfdz 
(if  at  reference  height  or  if  entrainment  is  not  considered)  and  dqdz  and/or 
dtdz. 

If  the  current  height  is  the  reference  height  or  entrainment  is  not  considered, 
(1/f)  df/dz  is  included  as  an  unknown  to  be  determined,  and  a  form  of  the 
Clausius-Clapeyron  equation  is  included  in  forming  the  matrix  equation.  If  the 
current  height  is  not  the  reference  height  and  entrainment  is  considered, 
(1/f)  df/dz  is  omitted  from  the  list  of  unknowns,  and  the  Clausius-Clapeyron 
equation  is  omitted.  The  variable  ntlpru  (Next  To  Last  Previous  Row  Used) 
is  defined  to  facilitate  using  much  of  the  same  code  for  both  cases. 

Unless  at  the  reference  height  or  in  a  situation  where  program  terminates 
abruptly  in  this  loop,  subroutine  drogesup  must  be  called  when  i  =  1  (because 
some  values  calculated  only  when  i  =  1  are  also  needed  in  subsequent  calls). 

7,4.10  Subroutine  iduncan 

Subroutine  iduncan  (that  has  formal  parameters  rlcm,  r2cm,  M2drops,  and 
TotDrops)  determines  the  number  of  droplets  per  cm^  having  radii  between 
rlcm  and  r2cm.  This  is  done  by  evaluating  a  closed-form  solution  to  a  sum 
of  definite  integrals  of  the  Duncan-Low  n(r)  density  function. 

The  formal  parameters  rlcm  and  r2cm  are  droplet  radii  in  centimeters.  The 
equivalent  radii  (rl  and  r2)  in  micrometers  are  used  when  passing  values  to 
function  xmeax. 

Formal  parameters  -  output  (real*16): 
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M2drops  -  concentration  of  mode  2  drops  (number  per  cm^)  having  radii 
between  rlcm  and  r2cm 

TotDrops  -  total  concentration  of  drops  (number  per  cm^)  having  radii 
between  rlcm  and  r2cm 

7,4A1  Subroutine  Illinois 

Subroutine  Illinois  uses  an  Illinois  iterative  process  (a  modified  Regula  Falsi 
method)  [9]  for  computing  a  solution  of  f(x)  =  fgiven. 

Formal  parameters: 

f  -  name  of  real*8  function,  which  must  be  available  as  a 

separate  routine  (input) 

fgiven  -  given  value  of  f(x)  (input;  reafS) 

X  -  root  of  f(x)  -  fgiven  =  0  (output;  real*8) 

xO  and  xl  -  two  approximations  of  x  such  that  f(xO)  -  fgiven  and 

f(xl)  -  fgiven  are  of  opposite  sign.  (These  values  may  be 
changed  by  routine.)  (input;  reaTS) 

delta  -  convergence  criterion:  If  change  in  value  of  estimate  of  root 

is  less  than  or  equal  to  delta,  then  root  is  considered  to  be 
found,  (input;  real*8) 

valid  -  validity  flag  (output;  integer): 

=  —  2  indicates  solution  was  not  found;  cause  unknown. 

=  —  1  indicates  error:  This  failure  occurs  when  the  input 
values  of  xO  and  xl  are  such  that  f(xO)  -  fgiven  and  f(xl)  - 
fgiven  have  the  same  sign  or  are  equal,  and  the  routine’s 
attempt  to  define  acceptable  initial  estimates  fails. 

=  1  indicates  valid  solution. 
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7.4.12  Subroutine  Inpcbloc 

The  cluster  sphere  radius  is  calculated  without  including  entrained  matter.  The 
resulting  value  is  stored  as  Pradius  in  common  block  /entrcomm/. 

If  the  initial  drop  size  distribution  was  a  bimodal  distribution,  the  mode  2  drop 
concentration  in  the  current  size  class  is  calculated  (if  necessary)  and  written 
into  a  binary  file  for  possible  use  by  another  program.  This  additional 
information  can  be  used  for  further  treatment  of  data  if  the  complex  indices  of 
refraction  for  drops  belonging  to  mode  2  are  assumed  to  be  different  from 
those  for  drops  belonging  to  mode  1 . 

If  a  drop  radius  decreases,  it  is  not  allowed  to  decrease  by  more  than  half  of 
the  difference  between  the  current  drop  radius  and  the  dry  radius  for  this  size 
class. 

7.4.13  Subroutine  interg 

Subroutine  interg  (that  has  formal  parameters  xl,  x2,  xi,  yl,  y2,  and  yi)  uses 
geometric  interpolation  to  determine  the  value  of  yi  for  the  point  (xi,yi)  when 
xi  and  the  points  (xl,yl)  and  (x2,y2)  are  given. 

7.4.14  Subroutine  inter! 

Subroutine  inter2  (that  has  formal  parameters  xl,  x2,  xi,  yl,  y2,  and  yi)  uses 
linear  interpolation  after  transformation  of  x  values  to  determine  the  value  of 
yi  for  the  point  (xi,yi)  when  xi  and  the  points  (xl,yl)  and  (x2,y2)  are  given. 
X  values  are  transformed  to  s  =  ln(l  -  x^)  before  linear  interpolation  is 
performed.  Y  values  are  not  transformed. 

7.4.15  Subroutine  inter! 

Subroutine  interS  (that  has  formal  parameters  xl,  x2,  xi,  yl,  y2,  ya,  yd,  and 
yi)  uses  linear  interpolation  of  transformed  values  to  determine  the  value  of  yi 
for  the  point  (xi,yi)  when  ya,  yd,  xi,  and  the  points  (xl,yl)  and  (x2,y2)  are 
given.  X  values  are  transformed  to  s  =  ln(l  -  x^)  before  linear  interpolation 
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is  perfonned.  Y  values  are  transformed  to  t  =  [(y  —  ya)/yd]‘-*  before  linear 
interpolation  is  performed.  After  linear  interpolation  is  performed  in  (s,t) 
space,  the  interpolated  value  for  t  is  subjected  to  an  inverse  transformation  to 
obtain  yi. 

7.4.16  Function  Inexpval 

The  real*8  function  Inexpval  (that  has  formal  parameters  m,  rl,  and  r2)  returns 
the  value  of  the  definite  integral  from  rl  to  r2  of  r*  n(r)  dr  divided  by  the 
expected  number  of  drops  in  this  interval,  where  n(r)  is  the  bimodal  lognormal 
density  function  with  parameters  given  in  common  block  /lognor/. 

7.4.17  Subroutine  Inndrops 

Subroutine  Itmdrops  (that  has  formal  parameters  rl,  r2,  M2drops,  and 
TotDrops)  returns  the  value  of  the  definite  integral  from  rl  to  r2  of  n(r)  dr, 
where  n(r)  is  the  bimodal  lognormal  density  function  with  parameters  given  in 
common  block  Aognor/.  Because  N1  and  N2  are  included  in  appropriate 
calculations,  the  returned  value  for  TotDrops  is  the  expected  total  number 
concentration  of  droplets  in  the  interval  from  radius  »  rl  to  r2,  and  the 
returned  value  for  M2drops  is  the  expected  part  of  that  concentration  attributed 
to  droplets  in  mode  2  of  the  bimodal  distribution. 

7.4.18  Function  LognMomt 

The  real*8  function  LognMomt  (that  has  formal  parameters  m,  rmode,  sigma 
a  in  equation  (1),  rlow,  and  rhigh)  returns  the  value  of  the  definite  integral 
from  rlow  to  rhigh  of  r*  n(r)  dr,  where  n(r)  is  the  lognormal  density  function 
defined  as 


^  exp  [-(Log  r  -  Log  rmodef  /  (2o^)] 
(Ln  10)  r  a 


(1) 


Ln  is  the  natural  logarithm  function,  and  Log  is  the  base  10  logarithm  function. 
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7.4,19  Function  LWCcUp 


7.4.20 


7.4.21 


7.4.22 


7.4.23 


The  real's  function  LWCcUp  (that  has  formal  parameters  nodmode,  nsizes,  pi, 
and  liioacgs)  returns  the  volume  of  condensed  matter  (cmVm^)  contained  in 
dr(q)s  having  radii  between  0.25  and  23.5  /tm. 

Function  LWCdiscr 

The  real's  function  LWCdiscr  (that  has  formal  parameters  nodmode,  nsizes, 
pi,  riioacgs,  and  writ)  returns  the  volume  of  condensed  matter  (cmVcm^) 
contained  in  aU  of  the  modeled  drops. 

Subroutine  LWCgrad 

Subroutine  LWCgrad  calculates  vertical  gradients  of  total  drop  volume  (Uquid 
water  content)  for  all  drops  and  for  only  those  drops  having  radii  between  0.25 
and  23.5  /tm,  inclusive. 

Subroutine  ndropi 

Subroutine  ndropi  stores  initial  numbers  of  droplets  in  an  array  that  wUl  not  be 
changed  after  initialization. 

Subroutine  PADLIMES 

During  a  given  iteration  for  a  given  height  level,  a  system  of  Unear  equations 
that  can  be  expressed  as  a  matrix  equation,  is  generated.  Subroutine 
PADLIMES  (figure  1)  solves  the  matrix  equation.  To  increase  efficiency, 
subroutine  PADLIMES  was  tailored  to  take  advantage  of  the  characteristics  of 
the  matrix  equations  this  subroutine  is  required  to  solve:  a  large  part  of  the 
coefficient  matrix  in  each  of  these  equations  is  an  identity  matrix  composed  of 
the  upper  left-hand  nsizes  X  nsizes  elements  of  the  whole  coefficient  matrix. 
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7.4.24  Subroutine  partiDSD 

Subroutine  partiDSD  partitions  a  continuous  drop  size  distribution  to  define  a 
discrete  distribution.  The  discrete  distribution  is  defined  to  have  a  finite 
number  of  size  classes  that  are  defined  for  sections  that  have  default  boundaries 
or  boundaries  set  by  user.  A  representative  radius  (mean  volume  radius)  is 
calculated  for  each  section.  (This  subroutine  is  used  if  sizmode  is  4  or  5.) 

7.4.25  Function  pnormal 

The  real*8  function  pnormal  (that  has  one  formal  parameter,  t)  calculates  the 
probability  that  a  normally  distributed  variable  with  mean  0  and  standard 
deviation  1  will  have  a  value  between  negative  infinity  and  t. 

7.4.26  Function  psatf 

The  real*8  function  psatf  (that  has  one  formal  parameter,  tzero)  uses  a  function 
from  Lowe  and  Ficke  [10]  to  calculate  the  saturation  vapor  pressure  over  water 
when  tzero,  the  temperature  (K),  is  given. 

The  range  of  validity  for  this  function  is  from  -50  to  +50  ®C. 
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Mubroutina  FAOLIMES  (cr,  cb,  y,  x,  ircr«  ircb,  nt,  nunk) 


Partially  Diagonalizad  Llnaar  Matrix  Equation  Solvar 
Solvaa  linaar  matrix  aquation  c  x  >  y 

whara  tha  qivan  coaffioiant  matrix  c  is  such  that  tha  uppar 
laft-hand  nt  x  nt  sub-matrix  is  alraady  an  idantity  matrix. 


* 
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1 

* 

I 

nt  X  nt  : 

cr 

* 

1 

idantity  : 

matrix. 

* 

1 

matrix  : 

dim 
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1 
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(nt,nb) 
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cb  matrix, 
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dimensioned  (nb,nunk) 

* 
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—  1 
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_  _  I 


Thara  is  no  naad  to  waata  computar  mamory  by  allocating  memory 
for  tha  antira  c  matrix.  Znstaad,  two  of  its  sub-matrices, 
or  and  cb,  arc  used. 

Formal  paramatarsi 

er  -  nt  by  nb  array  of  coafficiants  on  left  hand  side 
(sea  above)  (raal*8;  input;  will  ba  destroyed) 
cb  -  nb  by  nunk  array  of  coafficiants  on  Isft  hand  side 
(saa  abova)  ( realms ;  input;  will  ba  destroyed) 
y  -  vector  (column  matrix;  length  nunk)  of  constants  on 
right  hand  side  (real*8;  input;  will  ba  destroyed) 

X  -  vector  (column  matrix;  length  nunk)  to  contain 
solution  (real*8;  output) 

ircr  -  actual  cow  dimension  of  cr  array  as  dimensioned  by 
calling  program  (integer;  input) 
ircb  -  actual  row  dimension  of  cb  array  as  dimensioned  by 
calling  program  (integer;  input) 
nt  -  order  of  tha  identity  matrix  located  in  upper  left-hand 
corner  of  c  matrix  (see  above)  (integer;  input) 
nunk  -  number  of  unknowns  (integer;  input) 

Definition: 

nb  -  number  of  rows  (or  columns)  not  occupied  at  all  by 
the  previously-mentioned  (nt  x  nt)  idantity  matrix 
(see  above)  (integer;  calculated  as  nb  >  nunk  -  nt) 

Pointer  arrays: 

pointr  -  Subscript  is  row;  value  is  column, 
pointc  -  Subscript  is  column;  value  is  row. 
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A******************** *4* A********************* **************** ********** 


Figure  1.  Descriptive  comments  of  subroutine  PADLEVCES. 
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7.4,27  Function  rl6toi 


The  real*16  function  rl6toi  (that  has  formal  parameters  base  and  power)  raises 
a  real*16  number  to  an  integer  power  if  the  integer  is  small.  Repeated 
multiplications  are  used  in  a  way  that  preserves  mathematically  expected 
behavior  of  a  real  number  (base)  raised  to  an  integer  power  (exponent): 

If  the  base  is  positive 

The  function  returns  a  positive  value. 

If  the  base  is  negative 

And  the  exponent  is  even  (including  zero),  the  function  returns  a  positive 
value. 

And  the  exponent  is  odd,  the  function  returns  a  negative  value. 

If  the  base  is  zero 

And  the  exponent  is  positive,  the  function  returns  zero. 

And  the  exponent  is  negative  or  zero,  the  calling  program  is  aborted  with 
an  error  message. 

Formal  parameters 
base  -  base  (real*16) 
power  -  exponent  (integer) 

7.4.28  Function  rhf 

The  real*8  function  rhf  (containing  one  formal  parameter,  radius)  returns 
relative  humidity  as  a  function  of  droplet  radius.  Subroutine  Illinois  calls  this 
function  in  the  process  of  solving  for  the  droplet  radius  when  relative  humidity 
is  given. 

7.4.29  Function  ritfrdry 

The  real's  function  rhfirdry  (containing  one  formal  parameter,  radry)  returns 
relative  humidity  as  a  function  of  dry  radius.  The  subroutine  Illinois  calls  this 
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function  in  the  process  of  solving  for  the  dry  radius  when  relative  humidity  is 
given. 


7»4.30  Subroutine  rhfunc 

Subroutine  rhfunc  returns  relative  humidity,  activity  of  water,  and  curvature 
correction,  given  the  radius  of  the  droplet  and  other  variables  in  common. 

7.4.31  Subroutine  setnab 

Subroutine  setnab  sets  values  of  the  parameters  alphal,  alpha2,  betal,  beta2, 
Nl,  and  N2  (stored  in  labeled  common)  according  to  Duncan  and  Low.  [11] 
The  values  used  are  appropriate  when  radius  is  expressed  in  micrometers.  This 
subroutine  has  one  formal  parameter  visib,  which  is  the  visibility  in  kilometers 
(real*8).  The  value  must  be  supplied  by  the  calling  program  unit. 

7.4.32  Subroutine  setnrs 

Subroutine  setnrs  sets  parameters  for  a  bi  modal  lognormal  distribution 
according  to  a  rural,  maritime,  or  urban  aerosol  model  of  Shettle  and  Fenn.  [3] 

Subroutine  setnrs  requires  input  values  for  all  three  formal  parameters: 

fzero  -  relative  humidity  (1  =  saturated)  (real*8) 
visibility  -  meteorological  range  in  km  (real*8) 

SFmodel  -  the  Shettle-Fenn  model  used  (integer): 

1  -  rural 

2  -  maritime 

3  -  urban 

Subroutine  setnrs  uses  the  interpolation  routines  inter2,  inter3,  and  interg. 
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7.4,33 


Subroutine  slapsrat 


Subroutine  slapsrat  defines  a  coefficient  and  the  right-hand  side  for  cluster 
sphere  lapse  rate  equation.  (The  name  of  this  subroutine  comes  from  Sphere 
LAPSe  RATe.) 

In  part  1,  a  few  initial  values  are  stored. 

Part  2  calculates  quantities  in  an  equation  that  results  from  the  first  law  of 
thermodynamics,  the  ideal  gas  law,  a  hydrostatic  approximat'on  for  pressure, 
expressions  involving  entrainment,  an  expression  intended  to  model  (very 
i*oughly)  diffusion  involving  air  above  the  cloud,  and  a  linear  term  that  helps 
in  carve  fitting  and  might  be  considered  to  represent  other  effects.  This 
equation  is  to  be  solved  with  others  for  vertical  gradients  experienced  by  the 
cluster  sphere. 

7.4.34  Subroutine  solver 

Subroutine  solver  has  formal  parameters  fgiven  (real'8),  rwet  (real*8),  and 
sizeclass  (integer).  Subroutine  solver  calls  subroutine  Illinois  to  compute  the 
droplet  radius  (rwet)  at  which  the  relative  humidity  in  equilibrium  with  the 
droplet’s  surface  equals  the  input  value  of  fgiven.  The  formal  parameter 
sizeclass  is  the  size  class  of  the  droplet. 

7.4.35  Subroutine  solvrd 

Subroutine  solvrd  has  form  d  parameters  rwet  (rear8),  radry  (real*8),  and 
sizeclass  (integer).  Subroutine  solvrd  calls  subroutine  Illinois  to  compute  the 
dry  particle  radius  (radry)  for  which  the  relative  humidity  in  equilibrium  with 
the  surface  of  a  droplet  having  a  wet  radius  equal  to  the  input  value  of  rwet 
equals  t.ie  relative  humidity  value  of  fzero,  which  is  available  through  labeled 
common  and  is  used  as  an  input  value.  Again,  the  formal  parameter  sizeclass 
is  the  size  class  of  the  droplet.  The  variable  rdry  of  the  /params/  common 
block  is  set  equal  to  the  value  determined  for  the  formal  parameter  radry. 
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7.4.36  Subroutine  wrfiLs 

Subroutine  wriiles  reads  stored  data  from  unformatted  file  COMB.lile  and 
writes  the  same  data  into  individual  files,  most  of  which  are  in  plotit  format. 

7.4.37  Function  xmeax 

The  real*16  function  xmeax  (which  has  formal  parameters  areS,  m,  xlreS,  and 
x2re8)  returns  the  value  of  the  definite  integral  from  xlreS  to  x2re8  of  x"* 
exp(are8*x)  dx.  Are8,  xlre8,  and  x2re8  are  real*8.  m  is  an  integer.  The 
value  returned  by  this  function  is  real*16. 

7.5  Modifications  Needed  for  Use  on  a  Cray  System 

Global  differences  between  a  version  that  runs  on  an  HP  9000  system  (this 
version)  and  one  that  runs  on  a  Cray  system  should  include: 


HP  9000 

■Cray 

dabs 

abs 

dexp 

exp 

dfloat 

float 

dlog 

alog 

dsign 

sign 

dsqrt 

sqrt 

sngl 

-(omit)- 

The  changes  are  appropriate  because  REAL*8  precision,  used  for  most  of  the 
floating-point  variables,  is  considered  double  precision  when  using  Fortran  77 
on  an  HP  9000  computer  but  is  considered  single  precision  when  using 
Fortran  77  on  a  Cray  computer.  Floating-point  constants  are  specified  as 
single  precision  for  the  Cray  version,  such  as  by  replacing  d  (in  this  version) 
with  e  (for  Cray  version).  Floating-point  constants  in  this  version  are 
expressed  in  a  form  that  facilitates  global  substitution  to  change  from  double 
precision  to  single  precision. 
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Acronyms  and  Abbreviations 


ARL  Army  Research  Laboratory 

ASCII  American  Standard  Code  for  Information  Interchange 

cgs  centimeter-gram-second 

EECSwave  Extinction  Etc.  Coefficients  for  8  WAVElengths 

HP  Hewlett-Packard 

LN  natural  logarithm  of 

MACACASM  Microphysics  and  Cluster  Ascent  Cloud  and  Subcloud 
Model 

nm  nanometers 

rhs  right-hand  side 

sr  steradian 
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